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SECTION  1 
INTRODUCTION 


Currently  the  weapon  environment  codes  ROSCOE  and  WEPH  include  a 
model  for  the  dust  that  is  entrained  by  a near-surface  nuclear  fireball, 
which  was  developed  originally  for  the  RANG  code.  Since  the  develop- 
ment of  the  original  model,  additional  theory  and  new  experimental  data 
from  high  explosive  (HE)  nuclear  simulation  test  events  have  become 
available  on  dust.  This  report  presents  new  and  improved  models  for  dust 
and  propagation  effects  developed  from  the  new  information  for  the  WEPH 
code.  Although  developed  for  the  WEPH  code,  these  models  are  applicable 
to  other  engineering  codes  modeling  the  effects  of  dust  on  radar  and 
communications  systems. 

We  begin  by  discussing  possible  areas  for  new  and  improved  dust 
models.  We  then  choose  a number  of  these  areas  for  detailed  development 
in  this  report. 

The  Mie  theory  is  used  to  calculate  the  extinction  and  backscatter 
coefficients  for  the  dust  particles.  The  Mie  theory  gives  the  exact 
solutions  for  the  scattering  and  absorption  of  an  electromagnetic  wave 
incident  on  a uniform  spherical  particle.  The  solutions  are  given  as 
infinite  series  of  complex  terms.  The  current  Mie  program  in  the  systems 
codes  uses  a simple  forward  recursion  scheme  to  calculate  the  successive 
Mie  terms,  which  are  summed  until  convergence  is  reached.  However,  in 
some  cases  the  forward  recursion  technique  is  not  stable.  For  example, 
when  the  incident  radiation  wavelength  is  small  compared  to  the  dust 
particle  size,  a large  number  of  terms  are  required  in  the  Mie  scries, 
and  the  forward  recursion  scheme  can  lose  accuracy  at  the  higher  terms. 
For  the  largest  dust  particle  now  modeled — 10  cm — the  current  Mie  routine 
should  be  stable  up  to  about  10  to  20  GHz.  In  the  past  most  radar  and 


ff 
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communications  frequencies  were  lower  than  this,  but  some  current  sys- 
tems of  interest  have  frequencies  well  above  this  range. 

For  the  WOF  code,  a different  Mie  technique  was  developed  that  is 
stable  and  avoids  all  errors  inherent  in  the  forward  recursion  tech- 
nique. The  method  uses  a stable  backward  recursion  algorithm  and  uses 
a continued  fraction  scheme  to  accurately  evaluate  the  starting  high- 
order  term.  This  technique  is  valid  for  all  orders  and  can  thus  be  used 
for  any  frequency  and  any  particle  size.  The  current  Mie  routine  can 
be  replaced  by  the  new  improved  routine.  The  new  routine  is  compact; 
the  computer  running  time  is  only  slightly  longer  than  the  old  routine, 
but  does  require  about  300  to  400  extra  storage  locations. 

The  WFPII  code  calculates  only  the  attenuation  of  radiation  passing 
through  a dust  region.  Radar  codes  such  as  RANC  and  ROSCOE  also  cal- 
culate the  backscatter,  ie,  clutter,  due  to  the  dust  regions.  The  cur- 
rent model  formulation  for  the  clutter  power  due  to  scattering  from  a 
dust  region  is  the  following  : 


C = 


<4">  VeA 


G-r  (6  ,<t>)GR(9  ,<{>) 
4itR2 


0UN_.,  sin  4 dR  d6  d$ 

D 1 V 


where 

P 

F 


= peak  transmitter  power  (watts) 

= pulse  compression  ratio,  or  pulse  integration  improve- 
ment factor 

R = range  to  the  scattering  region  (m) 

9,<})  = angular  coordinates  with  respect  to  the  radar 

pointing  direction  (rad) 

G ( G , 4>)  = transmitter  antenna  gain  function  (including  system 
loss  factors) 

Gp(9>4>)  = receiver  antenna  gain  function  (including  system  loss 
factors) 


N 


PT 


= volumetric  number  density  of  dust  particles  (m  ' ) 

i 

= average  backscatter  cross  section  per  dust  particle  (m  ) . 


This  current  backscatter  formulation  implicitly  assumes  that  the  optical 
thickness  of  the  dust  cloud  is  very  thin,  the  total  one-way  attenuation 
being  less  than  about  1 to  2 dB.  Tor  typical  fireball  dust  clouds  and 
radar  frequencies,  the  optical  thinness  assumption  is  valid.  But  the 
dust  attenuation  rapidly  increases  as  the  frequency  increases,  and,  for 
some  high-frequency  radars,  the  dust  clouds  may  no  longer  be  optically 
thin.  The  full  formulation  for  a region  of  any  optical  thickness  is 


W/ 

<p  d R 


GT(e,<j>)GR(e,<}0 


°bntv 


e sin  <j>  dR  d';  d-; 


where 


optical  depth  at  range  R (due  to  all  sources  of  attenuation). 


If  only  the  dust  particles  contribute  to  the  attenuation,  then 


0HNTVdR' 


where 


= average  extinction  cross  section  per  dust  particle  (m  )• 


The  additional  term  e in  the  clutter  power  equation  accounts  for  the 
two-way  attenuation  of  the  radar  beam  to  and  from  the  scattering  volume. 

Although  the  extended  formulation  accounts  for  finite  optical  thick- 
nesses, it  is  still  a single  scatter  solution  (as  is  the  current  back- 
scatter calculation).  The  contributions  due  to  radiation  that  is  multiple 
scattered  are  ignored.  If  the  extinction  is  due  almost  entirely  to 
scattering  (little  or  no  absorption),  then  for  an  optically  thick  region 
the  backscatter  due  to  the  multiple  scattered  radiation  can  easily  dominate 
the  single  scattered  radiation.  To  first  order,  multiple  scattering  can  be 
accounted  for  by  adding  a build-up  factor,  B (0,<J>),  to  the  single  scatter 
integration  (see,  for  example,  Reference  1).  B^(0,<|>)  depends  upon  the 

optical  thickness  of  the  dust  cloud  in  the  direction  corresponding  to 
angular  coordinates  0 , cp , and  upon  the  dust  particle  sice  distribution. 


a 1 1 c nu.it  i on 


It  would  be  a relatively  easy  task  to  include  the  e 
term  in  the  backscatter  formulation.  It  would  require  considerably 
more  effort  to  include  the  multiple  scatter  build- tip  factor,  U 

The  current  dust  model  assumes  that  the  site  distribution  of  the 
dust  particles  can  be  described  by  a power  law  probability  distribution 
with  a power  exponent  of  4.  four  is  a typical  value  for  dust  particles 
from  loose  unconsolidated  soils  such  as  desert  alluvium,  bust  generated 
from  a nuclear  cratering  explosion  in  rock  and  cohesive  soils  has  a power 
exponent  of  about  3.5.  Very  fine  soils  may  have  a value  near  3. 

Another  very  common  site  distribution  is  the  log-normal  distribu- 
tion. This  distribution  is  often  used  for  small  particles  and  for  con- 
densates such  as  recondenscd  weapon  debris  and  water  droplets.  Since  the 
atmospheric  nuclear  test  ban,  a number  of  high  explosive  (HI  J surface 
or  shallow  buried  devices  have  been  detonated  to  simulate  nuclear  ex- 
plosions. Analysis  of  the  dust  particle  site  statistics  from  a number 
of  these  lit:  events  and  from  previous  nuclear  test  events  indicates  that 
the  smaller  dust  particles  iollow  a log-normal  d i st r i but  ion , while  the 
larger  particles  follow  a (tower  law  distribution.  A hybrid  site  dis- 
tribution can  be  defined  which  consists  of  a log-normal  joined  to  a power 
law  distribution. 

The  complex  index  of  refraction  determines  the  scattering  and  absorp- 
tion properties  of  a dust  particle.  The  present  dust  model  allows  for 
three  different  soil  types  (wet  clay,  dry  sand,  and  soil  with  an  ice 
coatingj . A built-in  table  gives  the  refraction  index  for  each  soil 
type.  These  indices  arc  taken  as  constants,  independent  of  the  incident 
radiation  frequency.  for  optically  thin  dust  regions  and  system  fre- 
quencies less  than  about  10  ('.lie,  a constant  index  of  refraction  is  a 
reasonable  approximation.  Tor  thick  regions  or  higher  frequencies, 
the  frequency  dependence  of  the  index  of  refraction  becomes  important. 

The  present  model  can  be  extended  to  include  the  frequency  dependence. 

I he  current  dust  model  considers  only  the  dust  entrained  in  the  fire- 
ball cloud.  In  addition  to  the  fireball  cloud,  there  are  two  other 
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regions  that  contain  dust — the  stem  and  pedestal  regions.  The  stem  is 
the  column  of  dust  reaching  from  the  ground  up  toward  the  main  cloud 
(the  "stem"  of  the  mushroom  cloud).  The  pedestal  is  that  relatively  low 
lying  dusty  region  extending  outward  along  the  surface  from  ground  zero. 
The  pedestal  is  formed  by  a combination  of  the  thermal  pulse  from  the 
fireball  and  the  outrunning  shock  wave.  There  are  photographic  data  from 
nuclear  test  events  on  the  growth  histories  of  the  stem  and  pedestal  re- 
gions. In  addition  there  have  been  many  theoretical  investigations  of 
both  regions.  There  is  now  sufficient  information  available  for  the 
development  of  a dust  model  for  these  two  subsidiary  dust  regions. 

Consider  a radiation  ray  path  through  a dust  region.  If  the  dust  is 
distributed  nonuniformly  in  the  region,  then  the  attenuation  and  back- 
scatter  effects  will  be  accompanied  by  fluctuations  about  the  mean  values, 
ie,  scintillations.  For  instance,  at  10  GHz  in  the  U11F/SHF  transmission 
experiment  on  the  dust  cloud  from  the  I.HCF  THROW  HH  explosion,  fluctua- 
tions of  ±8  dB  about  the  mean  attenuation  of  20  JB  were  observed  at  early 
times.  There  were  corresponding  fluctuations  in  the  phase.  The  HIGH 
THROW  cloud  was  highly  turbulent  during  these  observed  fluctuations. 

There  has  been  a large  theoretical  effort  devoted  to  the  transport  of 
radiation  through  regions  of  randomly  varying  dielectric  constant.  This 
corresponds  to  transport  through  turbulent  air.  In  the  field  of  nuclear 
effects,  transport  through  striated  plasma  regions  has  been  studied  and 
modeled.  However,  transport  through  a turbulent  dust  region  is  signif- 
icantly more  complicated.  Very  little  theoretical  work  exists  on  the 
distribution  of  dust  in  a turbulent  region.  Dust  particles  both  absorb 
and  scatter  the  incident  radiation;  the  scattering  is  typically  wide- 
angle  scattering. 

The  level  of  effort  required  to  develop  dust  scintillation  models  is 
not  known,  but  may  be  considerable.  A theoretical  investigation  would  re- 
quire the  mass  loading  and  sizes  of  the  dust  inhomogenit ies , and  their 
number,  location,  and  velocity  probability  distributions.  All  these 
parameters  must  be  known  as  a function  of  time.  Then  the  transport 
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through  a random  wide  angle  scattering  region  has  to  be  calculated. 

Rather  than  carrying  out  a full  theoretical  study,  it  may  be  possible  to 
develop  crude  first-order  models  based  on  experimental  data  and  simplified 
transport  physics. 

In  this  report  we  choose  the  following  four  areas  for  detailed  de- 
velopment : 

• Mie  scattering  theory 

• Size  distributions  of  dust  particles 

• Complex  index  of  refraction  of  dust  particles 

• Dust  models  for  the  stem  and  pedestal  regions. 

The  backscatter  and  scintillation  areas  are  not  discussed  further  in  this 
report . 


SECTION  2 

MIE  SCATTERING  THEORY 


The  Mie  theory  gives  the  exact  solutions  for  the  scattering  and  ab- 
sorption of  an  electromagnetic  wave  incident  on  a uniform  spherical 
particle.  The  solutions  are  given  as  infinite  series  of  complex  terms. 
Many  techniques  have  been  developed  to  numerically  evaluate  the  Mie  equa- 
tions. The  simplest  technique  is  to  use  a forward  recursion  scheme  to 
evaluate  the  successive  terms.  This  is  the  technique  used  in  the  present 
KHP11  code.  However,  in  some  cases  the  forward  recursion  technique  is 
not  stable.  To  avoid  the  instability  of  the  forward  recursion  tech- 
nique, many  backward  recursion  schemes  have  been  developed.  For  the  WOE 
code,  Reference  1,  a new  backward  recursion  Mie  code  was  developed.  This 
Mie  code  is  stable,  accurate  at  all  orders,  and  can  be  used  for  any  fre- 
quency and  particle  size.  The  following  description  of  the  new  Mie  cal- 
culation is  largely  taken  from  the  WOE  code  documentation.  Reference  1. 
Appendix  A is  a listing  of  the  new  Mie  computer  routines. 

We  first  present  the  formulas  for  the  Mie  solution  (see  any  standard 
text  for  a derivation)  and  then  the  method  used  to  solve  the  equations. 
Define  the  following  quantities: 

a = diameter  of  the  spherical  particle  (cm) 

A = wavelength  of  the  incident  radiation  (m) 

in-2  713 

a = 10  T 

-2 

= dimensionless  size  parameter  (the  factor  10  converts 
a from  cm  to  m) 


m = 


im 


1 


= complex  index  of  refraction  of  the  sphere  (note  that  here 


we  are  using  in  instead  of  n,  since  by  custom  n is  used  as 
the  order  in  the  Mie  formulas) 


| 
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= mot 


°SCA 


ABS 


BKS 


<SCA 


= scattering  cross  section  of  the  sphere  ( m“) 

2 

= absorption  cross  section  of  the  sphere  (nrj 

2 

= backscatter  cross  section  of  the  sphere  (m  ) 


= 10 


JSCA 


TT(a/2)' 


= scattering  efficiency  (ratio  of  scattering  cross  section  to 

4 

physical  cross  section)  of  the  sphere.  The  factor  10  converts 

2 2 

the  sphere  area  in  cm  to  m . 


<ABS 


= 10 


4 °ABS 


<EXT 


QBKS 


tt  (a/ 2) 

absorption  efficiency  of  the  sphere 

Q^BS  + Qsca  = extinction  efficiency 
. Or 


= 10 


BKS 


Tt(a/2) 

backscatter  efficiency 


-] 


S(0)  = scattering  function  (m“  sr  J).  Sic)  dl.  is  the  fraction  of  the 

incident  unpolarized  energy  per  unit  area  that  is  scattered 

into  solid  angle  dfi  centered  about  the  direction  that  makes 

an  angle  6 with  the  direction  of  the  incident  radiation 

(0  is  the  scattering  angle).  Note  that  o „ = 4ttS(tt). 

dKo 


Currently  the  WLPH  code  only  uses  the  extinction  parameters  (to 
calculate  total  attenuation).  We  include  the  scattering  and  absorp- 
tion parameters  and  the  scattering  function  in  our  present  discussion 
of  the  Mie  solutions  for  completeness,  and  to  allow  for  use  in  other 
DNA  engineering  codes. 
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I'he  equations  for  the  Mie  solution  are: 


‘'SCA  = 4 £ <2n  * fl  a„  I ^ * lb„|2] 

a n=l  L J 


<*txr  * T s <2"  ‘ » Re  (»„  * V 

a n=l 

(where  Re  signifies  the  real  part  of) 


^ABS  ^LXT  ' ^SCA 


oo 

= A i Z (n  + T)(-l)n  (i 


a ni=l 


a - b ) I4 
n n 1 


where 


SO)  = j { !s1  ce) |2  ♦ |s2(9)|2|  , 


a'T(Y)4'  (a)  - YY  (a)4(  (Y) 
n n n n 

Jn  = (Y)^(a)  - Y£n(a)Yn(Y) 


Yyn(Y)4'n(a)  - “^WyY) 

bn  = YY^Yj^Ca)  - ai;n(a)Yn(Y) 


00 

SI(9)  ' £ Win  I i)  {Vn(cos  6)  * Vn(cos  0)} 

s2(e>  * £ wfirTTy  {v„<cos  9>  * Vn<cos  9)} 

/ \1/2 

VZ>  = (t)  Jn*l/2<*> 

^Z>  =(t)  [Jn*l/2(Z)  + i(-])n  J-n-l/2^J 


(12) 


71  (COS  Uj  = I’  (COS  0) 

n n 


2 d 

(cos  0)  = COS  0 71  (cos  0)  - sin”0  jr  71  (cos  0)  . (13) 

i n J d cos  On 


The  J's  arc  spherical  Bessel  functions  of  complex  argument  and  half- 
integer order.  The  l”s  arc  Legendre  polynomials.  T and  f are  Riccati- 
Bessel  functions,  and  it  and  i are  associated  Legendre  polynomials.  De- 
fine an  arbitrarily  oriented  plane  containing  the  scattering  sphere  and 
the  incident  radiation.  Then  for  scattered  radiation  within  the  plane, 
the  complex  amplitude  function  Sj (0)  describes  the  scattering  for  an  in- 
cident plane  wave  with  vertical  polarization  (L  perpendicular  to  the 
scattering  plane;  S->(0)  is  for  horizontal  polarization  (E  parallel  to 
the  scattering  plane). 


As  might  be  expected  from  the  complexity  of  the  Mie  equations,  the 
numerical  evaluation  of  Q Q^xi’  ant*  S(Gj  for  a given  m and  a is  not, 
in  general,  a trivial  task.  The  terms  in  the  infinite  series  have  to 
be  evaluated  and  summed.  The  number  of  terms  that  have  to  be  evaluated 
before  the  series  converge  depends  primarily  upon  the  size  parameter  a. 
Roughly,  the  number  of  Mie  terms  required  is  1.5a;  for  large  particles 
and  small  wavelengths,  several  hundred  terms  are  often  required  before 
convergence. 


The  evaluation  of  the  various  orders  of  C , v , and  x are  straight- 

n n n 

forward.  We  can  use  the  well  known  recurrence  relations  of  Bessel  func- 
tions and  Legendre  polynomials  to  obtain 


L 

n 


(a) 


cn 


(a) 


,la) 


(14) 


with  initial  values 

f j(a)  = sin  a + i cos  a 

( j(a)  = cos  a - i sin  a 

TT  (COS  0)  = — r-  COS  9 7T  (COS  9)  - p 71  - (COS  0) 

n ' n - 1 n-1  n - 1 n-2 


(15) 

(16) 
(17) 
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•1 


\ ; 

t 


T (cos  0)  = COS  0 Tl  (cos  0)  - n _(cos  0) 
n L n n-2  J 


- (2n  - 1)  sin“0  n^^fcos  6)  + in  2(cos  G) 


(18) 


The  initial  values  are 


ti  (cos  0)  = 0 
o J 


i (cos  0)  = 0 
o 


1tj(COS  0)  = 1 


T j (COS  0)  = COS  0 


TF,(COS  0)  = 3 COS  u 

T-,(cos  0)  = 3 cos  (26) 

With  the  initial  values,  we  can  use  the  forward  recurrence  relations  to 
generate  the  required  terms  to  any  order.  The  forward  recursion  tech- 
nique for  these  three  functions  is  stable  and  accurate. 

To  complete  our  numerical  evaluation,  we  define  the  complex  function 


H1  ( Y \ 

VY>  = w 


(19) 


with  this  definition,  the  Mie  formulas  for  a^  and  can  be  written: 


(—  * ") 


Re  \ (u)j-  - Re 

In  ) 


{W“>! 


(20) 


/A  00  \ 

— + -k  (a)  - C . (a) 

\ m a/  n n-1  ' 


K<V>  • |)lte{t„(a)}  - (4 


ImA  (Y)  + -k  (a)  - t)  . 
\ nv  ' a)  n n-1 


(a) 


(21) 


The  primary  difficulty  in  the  evaluation  of  the  Mie  formulas  lies  in 


the  evaluation  of  An(Y).  Using  the  properties  of  the  Bessel  functions. 
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I 


T_ 


•• 


I 


i 


we  can  write  A (Y)  as 
n 


An0') 


n 

Y 


Jn-l/2 

Jn+l/2 


(Y) 

00 


(22.) 


Thus,  if  we  can  evaluate  the  Bessel  functions,  say  by  forward  recursion, 
A^ (Y j , can  be  evaluated.  Alternately,  we  can  use  the  recurrence  rela- 
tions for  the  ratios  of  the  Bessel  functions  and  write  the  recursion 
equation  for  A^(Y)  itself  as 


vv)  ■ -?*  (?-  vH'1 


with  initial  condition 


Ao(Y) 


cos  Y 
sin  Y 


(24) 


The  forward  recursion  technique  for  the  evaluation  of  A (Y)  is  very 
susceptible  to  error  in  at  least  four  cases  (Reference  2): 

• When  the  argument  is  small 

• When  the  argument  is  large,  requiring  a large  number 
of  orders 

• When  the  imaginary  value  is  larger  than  the  real  value 

• For  certain  anomalous  values. 

The  use  of  forward  recursions  to  generate  the  consecutive  orders  of 
Bessel  functions  is  a classic  example  of  unstable  numerical  methods. 

The  current  WhPlI  Mie  calculation  uses  an  asymptotic  analytic  formula 
whenever  the  argument  is  small  and  thus  avoids  the  first  error  case;  the 
others  still  remain,  however. 

Many  other  techniques  have  been  devised  to  generate  the  required 
Bessel  functions  or  ratios.  Most  techniques  involve  some  type  of  back- 
ward recursion.  The  values  of  the  Bessel  functions  or  the  ratios  arc 
evaluated  at  a high  order,  and  the  backward  recursion  relation  is  used 
to  evaluate  the  lower  orders.  The  backward  recursion  does  not  suffer 


rife 
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from  the  instability  of  the  forward  method.  However,  care  must  be  taken 
to  preserve  accuracy;  some  techniques  lose  accuracy  even  when  using 
double  precision  arithmetic.  Recently  Lentz  (Reference  2)  has  developed 
an  algorithm  for  evaluating  the  Bessel  functions  and  latios  that  elim- 
inates the  weaknesses  of  the  earlier  methods.  Lentz's  algorithm  uses  a 
new  technique  of  evaluating  continued  fractions  that  starts  at  the  be- 
ginning rather  than  the  tail  and  has  a built-in  error  check.  Using  the 
method,  any  A (. V _)  can  he  generated  completely  independently  of  all  pre- 
ceding values  with  high  accuracy.  Readers  are  referred  to  Lentz's 
article  for  details. 

We  use  Lentz's  method  to  generate  A (Y)  for  n of  order  « 1.5a  and 
then  use  the  backward  recursion  relationship. 


to  generate  all  lower  orders.  Using  the  forward  recursion  relations  for 

the  other  functions,  the  a and  b are  calculated  and  the  infinite  series 

n n 

summed  until  convergence.  In  almost  all  cases,  convergence  is  reached 
before  reaching  the  highest  precomputed  order  of  A (Y).  Otherwise 
Lentz's  method  is  used  to  generate  any  additional  needed  terms. 

Utilizing  the  Lentz  algorithm,  we  have  written  a very  compact  com- 
puter routine  that  evaluates  the  exact  Mie  equations  for  Q Qggs> 

Q. ......  (and  thus  Q ),  and  S(b).  The  running  time  is  quite  reasonable 

for  an  exact  calculation.  For  a = 1.2,  three  orders  are  required  for 
convergence,  and  the  running  time  is  1 millisecond  on  a GDC  7600  com- 
puter. lor  a = 100,  103  orders  are  required  with  a running  time  of  25 
milliseconds.  For  inclusion  in  the  Will'll  code,  a simpler  routine  has 

been  assembled  which  calculates  only  Q.......  and  QD1/C;  this  routine  runs 

LA  l dno 

considerably  faster  than  the  full  calculation  case.  Listings  of  both 
the  full  Mie  calculation  routine  and  the  reduced  routine  are  given 
in  Appendix  A. 


SECTION  3 

SIZE  DISTRIBUTIONS  OF  DUST  PARTICLES 


T 


The  two  most  common  probability  distributions  used  to  describe  par- 
ticulate (whether  dust,  smoke,  haze,  fog,  rain,  or  debris)  size  statistics 
are  the  power  law  and  the  log-normal  distributions.  The  power  law  prob- 
ability distribution  is 


where 


P (al  _ CP  - 1J  a~p 

P n ~ C P - 1 ) _ a- (p-i) 


min 


a . < a < a 

nun  — — max 


(26) 


max 


a = particle  diameter  (cm) 

Pp (a)  da  = fraction  of  the  particles  with  diameters  between 
n and  a + da 

p = power  law  exponent 


min 


= minimum  particle  diameter  in  the  distribution  (cm) 


a = maximum  particle  diameter  in  the  distribution  (cm). 


The  log-normal  distribution  is 

_2 


PL<a>  = 


Xn  a/ a 


Jin  S 


yfl-n  a £n  S 


0 < a < a> 


(27) 


where 

a^  = mean  particle  diameter  (cm) 

S = standard  deviation  parameter. 

The  two  distributions  given  are  probability  distributions,  ic, 
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/ I>  (aj  da  = / I’ L ( a ) da  = 1 

•L  . *'o 

nun 


N.  = total  number  of  particles  in  the  power  law  distribution 
N l = total  number  of  particles  in  the  log-normal  distribution. 


y»)  = 

fL<“)  = 'Vp1"’ 

are  the  number  distributions  for  the  two  cases,  where 

f^a)  da  = number  of  particles  in  the  power  law  distribu- 
tion which  have  diameters  between  a and  a + da. 


M.|.  = total  mass  of  particulates  (g) 

p = bulk  density  of  the  particulate  material  (g  cm  *) , 


6 Va  >p 


where 


* n 


71  ,3, 

6 Va  L 


a P (a)  da 
P 


P = 4 


a-fP-l)  _ a-(p-l) 
1 min  max 


1 ra4-p . a4-p_ 

^4  - p I max  min  _ 


p t 4 
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<) 


a^P.  (a)  da  = a'1  e 
L m 


T I kn  SI- 


133  j 


Condensates,  such  as  water  droplets  or  the  particles  reformed  from 
the  weapon  vaporized  material,  are  generally  well  described  by  the  log- 
normal size  distribution.  Some  typical  values  used  in  previous  parti- 
culate models  are  a = 1.4  x 10  ^ cm  (0.14  pm) , S = 1.9  for  recondensed 
m -4 

weapon  material  particulates  and  a = 5 x 10  cm  (5  pm),  S = 2.0  for 
water  droplets  within  a nuclear  cloud. 

Oust  and  crater  ejecta  part iculates  are  generally  better  described 
by  a power  law  distribution,  at  least  for  particles  greater  than  a few 
hundredths  of  a centimeter  in  diameter,  experimental  data  from  a number 
of  Ilf  (high  explosive]  and  nuclear  tests.  References  5 and  4,  indicate 
that  typical  values  for  the  power  law  exponent  are  p ~ 3 . 5 for  rock  and 
cohesive  soils  (such  as  clay  or  shale)  and  p ~ 4 for  loose  unconsoli- 
dated soils  (such  as  desert  alluvium  or  sand),  fine  soils  may  have  p - 5. 

The  log-normal  distribution  is  well  behaved  mathematically  at  both 
limits  of  small  and  largo  particle  diameters.  The  power  law  distribu- 
tion is  not  well  behaved  at  either  limit.  As  a . goes  to  zero  both  the 

nun  - 

probability  distribution  P (a  . ) and  the  total  number  of  particles  .\_,r, 

1 • p nun  1 IP 

go  to  infinity.  As  a x goes  to  infinity,  the  probability  distribution 

l’(a  ) goes  to  a well  behaved  zero,  but  goes  to  zero  for  a finite  M... 
max'  11  1 

There  are  several  techniques  employed  to  eliminate  mathematical  dif- 
ficulties at  the  endpoints.  The  most  common  technique,  and  the  one  we 
Implicitly  assumed  when  we  defined  the  power  law,  is  simply  to  cut  the 
distribution  off  at  some  lower  and  upper  limits,  a . and  a . Another 
technique  is  to  assume  an  upper  limit,  a but  to  attach  a mathematic- 
ally well  behaved  tail  at  the  lower  end  of  the  distribution.  This  lower 
limit  tail  can  be  some  artificial  mathematical  expression  chosen  for 
computational  convenience,  or,  as  is  often  the  case,  can  be  a fit  to  the 
actual  experimental  data  in  the  small  particle  limit. 
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SAI,  Reference  4,  reviewed  the  experimental  data  on  dust  particle 
sice  distributions  from  nuclear  test  events.  They  then  analyzed  in  de- 
tail the  more  recent  (and  more  complete!  dust  size  data  from  111:  tests. 
Their  conclusion  was  that  the  small  size  particles  have  a log-normal 
size  distribution,  while  the  larger  particles  huv  a power  law  size 
distribution.  The  size  division  between  small  and  large  particles  oc- 
curs at  a particle  diameter  of  about  0.018  cm  [180  pm).  Ke  will  develop 
the  equations  for  the  hybrid  size  distribution  of  a log-normal  small - 
limit  tail  attached  to  a power  law  distribution. 


The  hybrid  probability  distribution  is  given  by 
jLl'  I.' 


(C  P (a)  0 < a < a . 


PhW  = lc 


[34) 


C2Pp(a) 


a . < a < a 

nun  — — max 


where  and  C0  are  normalization  constants  which  ensure  that 


. a 


max 

P.,[a)  da  = 1 


(55) 


<3,1’,  (a  • ) 

1 L nun 


= C\P  [a  . ) 

> xx  v mi  tv' 


min' 


(36) 


ie,  the  constants  ensure  that  the  total  probability  distribution  is 

properly  normalized  to  unity  and  that  the  log-normal  and  the  power  law 

distributions  join  at  a . . C,  and  (h  arc  given  by 

J min  1 2 b '■ 


P (a  . ) 
p nun 

P.  (a  . ) 
L innr 


(57) 


C2 


P 

J1 


(a 

V 


min'* 

min^. 


-1 


where 


(58) 


I 

I 
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The  fractions  of  the  total  dust  mass  which  lie  in  the  two  segments  are 


lor  cohesive  soils  Reference  1 recommends  the  following  values  for 
tile  hybrid  distribution  parameters: 

a = D.tHU  cm  I JO  (.ml 

m 

S = J 

a = 0.0 IS  cm  (iso  urn) 

nun 

a = 1 cm 

max 

p =5.5  - 

The  current  IVHI'll  model  assumes  that  for  nuclear  bursts  the  size 
distribution  is  a power  law  distribution  best  represented  as  that  of 
unconsolidated  soils  and  takes 

i 

P (a)  = =-  a 4 , ie,  p = I | 

P a-'  - a ’ 

‘ min  max  j 

i 

a ■ =0.001  cm  (10  |.m|  i 

nun  ' j 

a =10  cm  . ] 

max 

The  10-cm  a value  is  for  surface  and  very  near  surface  (ie,  cratering i 

max  • j 

bursts.  As  the  burst  height  is  increased,  cratering  ceases,  the  lofting  ; 

power  of  the  nuclear  induced  winds  decreases » and  a is  also  assumed  to  j 

1 max 

decrease.  1 

he  can  improve  the  size  distribution  model  by  adopting  cither  \ 

a general  power  law  distribution  (arbitrary  p,  not  just  p = 4) , or  j 

the  general  hybrid  distribution.  As  we  shall  show  later,  the  attenua-  j 

tion  and  backscattcr  effects  of  the  hybrid  distribution  are  not  signi-  \ 

ficanliy  different  from  t'ne  corresponding  power  law  distribution  except  ; 

for  extremely  high  frequencies  (greater  than  about  100  (lllz).  for  the  i 

present,  we  adopt  the  generalized  power  law  distribution.  If  in  the  j 

future  the  higher  frequencies  become  important,  the  hybrid  distribu- 
tion can  be  easily  implemented  into  the  hiiPII  code.  The  changes  required 
in  the  code  to  generalize  the  fixed  p = 4 size  model  to  a general  p 
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are  trivial.  Appendix  B gives  a list  of'  the  generalized  model  equations 
and  a listing  of  the  generalized  I’HKOUP  subroutine,  which  calculates  the 
extinction  and  backscatter  cross  sections.  To  illustrate  the  variation 
of  effects  with  changing  size  distribution  parameter,  figure  1 shows  the 
, extinction  cross  section  per  gram  of  material  as  a function  of  p for  dry 

l 

f sandy  soil  tor  a number  of  frequencies. 

jj.  For  a given  amount  of  dust  mass,  what  is  the  difference  in  extinction 

I and  backscatter  properties  for  the  hybrid  and  the  corresponding  power  law 

j;  : size  distributions?  To  illustrate  the  differences  numerically,  we  adopt 

j the  following  parameter  values: 


Power  law  size  distribution  Hybrid  size  distribution 


a . 
nun 

= 0.001  cm  {10  pm  J 

a . = 0.018  cm  (ISO  pm J 

m i n 

a 

= 10  cm 

a - 10  cm 

max 

mu  x 

P 

= 3.5  and  4 . 0 fie, 

a = 0.002  (20  pm) 

m 

hard  rock  and  soil 

S =2 

part iclesj 

p =5.5  and  4.0 

For  both  distributions,  we  take 

Mr  = HZ'  g 

= --5  g cm"'1 

The  probability  distributions  are 

(7.906  x 10~8  a"'5*5 
pp(a)  = ! 

1 / -9  -4 

<3.000  x io  a 

The  hybrid  distribution  is 


p = 3.5 


0.001  < a < 10 


p = 4 . 0 


(46) 


PH(a)  = 


0.5750 


1 [*n  Ool" 

2 P.n  2 


1.643  x 10'7a'J‘5 


0 £ a 5 0.018 
0.018  < a < 10 


0.5751 


on  

1 0,002 
2 [ £n  2 


-8  -4 

2.204  x 10  a 


0 < a < 0.018 


0.018  < a < 10 


p = 3.5 


p = 4.0 


Figure  2 shows  the  power  law  and  hybrid  distributions.  In  addition 
the  entire  log-normal  probability  distribution  is  shown.  We  see  that 
the  log-normal  and  hybrid  distributions  have  many  particles  below  the 
small-limit  cutoff  of  the  power  law  distribution.  Both  the  power  law 
and  hybrid  distributions  have  more  large  particles  than  the  log-normal 
distribution. 

The  normalization  constants  for  the  hybrid  distribution  are 
( 0.9991  p = 3.5 

C = (48) 

| 0.9992  p = 4.0 


j 1.512  x 10' 
| 1.260  x 10' 


p = 3. 5 
p = 4.0 


The  (a  ) values  for  the  two  distributions  are 


j 4.950  x 10 
| 3.482  x 10 


-7 

3 

i 

cm 

p = 3.5 

f8 

3 

p = 4.0 

cm 
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Figure  2.  Dust  particle  size  probability  distributions. 
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/ 1 . 055  x 10"6  cm3  p = 3.5 

(a3)“  = | -73 

'1.992  x 10  cm  p = 4.0 


The  total  number  of  particles  in  10  grams  of  mass  is 


(N't)-  = i 


1.543  x 10  p = 3.5 


T P | 13 

V2.665  x 10 


2.665  x 10  p = 4.0 


(7.241  x 10  p = 3.5 


OM..  = 


"'I 


3.835  x io  p = 4.0 


The  number  distributions  for  the  two  cases  are 


f (a)  = (N,r)  P (a)  = 
u ^ ' I p p 


fl.220  x 104 5a"J'5 


I 4 -4 

'8. 295  * 10  a 


p = 3.5 

0.001  < a < 10 
p = 4 . 0 


4 . i(. ; • in 


) n ■ ■ 


(i  ■ ;i  • (>.  •*»  j ?• 


/ 

) 

{ 1 . 1 ill)  ■ 1C  a ' 


P : 3.1 


o.nui  i ■ io 


y:.)  - (V.iV  • 


C.200  • Id 


4 -4 

8.452  - 10  a 


r - <*>  i 

1 ' " in',  on': 

' I J 


(I  < a < 0. 0!  8 


0.01*  •'  a • 10 


We  see  that  although  the  probability  distributions  and  the  total  number 
of  particles  arc  different  for  the  power  law  and  the  hybrid  distributions, 
the  number  of  particles  versus  size  is  virtually  identical  for  particles 


greater  than  0.018  cm  (the  joining  point  between  the  log-normal  segment 
and  the  power  law  segment  in  the  hybrid  distribution). 


r 


t i 


£_^a) 


fH(a) 


0.9752 

p = 3.5 

1.019 

p = 4.0 

a > 0.018 


(56) 


So  any  differences  in  the  extinction  and  backscatter  properties  are  due 
almost  entirely  to  the  smaller  particles  (those  with  diameters  less  than 
0.018  cm  (180  pm)).  The  fraction  of  the  total  number  of  particles  with 
diameters  less  than  0.018  cm  and  the  fraction  of  the  total  mass  carried 
by  these  smaller  particles  are 


( 0.9993 

p = 3.5 

*-Vp  " 

1 0.9998 

a 

p = 4.0 

£ 0.018 

(57) 

(Vh  = 

J 0.9991 

) 

p = 3.5 

1 a 

£ 0.018 

(58) 

1 0.9992 

p = 4 . 0 
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( 3.278  x 10 

p = 5.5 

<fN>p  = 

) 

1 0.3139 

p = 4.0 

a £ 0.018 

(59) 
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(5.680  x io 

p = 3 . 5 

<Vll  = 

) 

1 0.3004 

p = 4.0 

a £ 0.018 

(60) 

Hence,  more  than 

99.9  percent  of  the  particles  are  small  particles. 

, but 

they  constitute 

only  about  3 to 

6 percent  of  the 

total  mass  for  a hard 

1 


rock  size  distribution  (p  = 3.5)  and  about  50  percent  for  a soil  dis- 
tribution (p  = 4.0). 


The  scattering  and  absorption  cross  sections  for  a particle  of 
diameter  a are  given  by 


. ..4  / a ' ~ 2 

0.,,..  = 10  u T Qt.,, , in 
SCA  \ 2 XSLA 


°ABS  = 1(}4tT(!)Q  ABS  ^ 

where  Qc„.  and  Q. ....  are  the  scattering  and  absorption  efficiencies, 
oLA  Ado 

respectively.  In  the  long  wavelength  limit,  ic,  where 


10"2Tra  . 

a = « 1 


A = wavelength  of  the  incident  radiation 


qsca  11  a 


^ABS  “ U 


so  that 


°SCA  “ a 


°ABS  a 3 • 


In  the  short  wavelength  limit,  ic,  where 


a » 1 , 


then  both  Qer..  and  Q,nc  approach  constant  values  and 

oLA  Ado 


°SCA  “ a 
°ABS  “ ^ * 

So  in  the  long  wavelength  limit,  it  is  the  largest  particles  that  domi- 
nate the  total  scattering  cross  section  for  the  whole  sice  distribution. 
In  this  limit,  the  absorption  cross  section  is  seen  to  be  proportional 
to  the  amount  of  mass  carried  by  the  particles.  Most  often  in  the  long 
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wavelength  limit,  absorption  dominates  scattering  so  that  t lie  extinction 
is  also  determined  by  the  amount  of  mass.  In  the  short  wavelength  limit, 
both  absorption  and  scattering  are  dominated  by  particles  with  diameters 
about  a particular  size.  This  dominant  size  is  generally  near  t = 1. 

Hence  if  the  long  wavelength  limit  obtains  for  the  small  particles, 
the  extinction  and  backscatter  properties  of  the  power  law  and  the  hy- 
brid distributions  are  virtually  identical.  The  scatter  is  determined' 
by  the  larger  particles,  and  both  numerical  distributions  are  essentially 
the  same.  The  absorption  contribution  of  the  smaller  particles  is  de- 
termined by  their  total  mass,  which  as  we  have  seen  is  negligible  for 
hard  rock  (p  = 5.5)  and  almost  identical  (51  and  50  percent)  for  the  two 
distributions  for  soil  (p  = 4.0).  If  the  long  wavelength  limit  fails 
for  any  of  the  small  particles,  then  the  two  distributions  will  have  sig- 
nificantly different  extinction  and  backscatter  properties. 

The  frequency  corresponding  to  u = 1 for  a particle  of  diameter 
0.018  cm  is 

f = 550  GHz. 

Hence,  for  frequencies  below  about  100  GHz,  the  small  particles  can  he 
considered  to  be  in  the  long  wavelength  limit,  and  the  power  law  and  hy- 
brid distributions  produce  essentially  identical  propagation  effects. 

The  simpler  power  law  formulation  can  be  used  without  loss  of  accuracy. 

If  the  frequency  of  interest  rises  much  above  100  GHz,  then  the  hybrid 
distribution  formulation  should  be  implemented. 


SECTION  4 


COMPLEX  INDEX  OF  REFRACTION 
OF  DUST  PARTICLES 


The  complex  index  of  refraction  determines  the  scattering  and  absorp- 
tion properties  of  a dust  particle.  In  the  current  Wf  PH  model,  there  are 
three  typical  soil  types  built  in  the  program  with  the  indices  given  in 
Table  1.  These  are  representative  values  for  soils  at  a frequency  of  3 
GHz  lS-bandJ . These  values  are  taken  as  constants,  independent  of  fre- 
quency. At  frequencies  much  lower  than  S-band,  the  dust  effects  are  gen- 
erally negligible.  It  was  also  assumed  that  the  primary  interest  was  for 
frequencies  below  about  10  GHz.  Therefore  a constant  index  of  refraction 
was  a reasonable  assumption.  But  with  the  addition  of  a stem  dust  model, 
the  dust  loading  may  be  high  enough  to  significantly  affect  frequencies 
lower  than  S-band.  In  addition,  radar  and  communication  systems  with 
frequencies  higher  than  10  GHz  are  becoming  more  common.  The  frequency 
range  of  interest  has  broadened  sufficiently  that  the  assumption  that  the 
index  of  refraction  is  independent  of  frequency  is  no  longer  valid. 


Table  1.  WEPH  soil  types. 


Soil  Type 

Index  of  Refraction 

Wet  clay 

3.5  - 0.4i 

Dry  sand 

1.5  - 0.025i 

Icp-covered  soil 

1.78  - 0.C024i 

Also,  allowing  only  three  choices  of  soil  types  is  unduly  restrictive. 
Ihcre  is  essentially  an  infinite  variety  of  soils,  with  a continuous  vari- 
ation of  the  index  of  retraction.  The  absorption  due  to  dust  particles 
is  proportional  to  the  imaginary  part  of  the  complex  index  of  refraction. 
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The  three  soil  choices  have  imaginary  indices  roughly  an  order  of  mag- 
' nitude  apart.  The  attenuation  for  a system  could  be  1,  10  or  '00  dB, 

depending  on  the  choice  made  for  the  soil  type.  A finer  d’  ion  of 
soil  types  (and  thus  attenuations)  is  called  for. 

Another  major  assumption,  implicitly  expressed,  is  that  the  soil 
index  of  refraction  can  also  he  used  for  the  index  of  refraction  l'or  the 
individual  dust  particles.  Most  soil  indices  of  refraction  are  derived  i 

from  measurements  made  on  soil  samples  which  have  been  removed  from  the  ’ 

field  and  carried  back  to  the  laboratory.  Soil  can  be  considered  to 
consist  of  three  components — solid  material,  water,  and  air.  The  soil 
indices  are  determined  by  the  relative  fractions  of  the  three  components, 
and  by  the  electrical  properties  of  the  solid  and  water  components.  In 
situ  soil,  laboratory  soil  samples,  and  individual  particles  will  all 
have  different  fractions  of  the  three  constitutive  components  and  will 
thus  have  different  indices  of  refraction.  There  are  various  mixing 
rules  available  which  relate  the  index  of  soil  in  one  state  to  the 

index  of  the  same  soil  in  another  state,  he  will  make  use  of  the  1 

somidisperse  model  of  Reference  5.  In  addition  to  relating  the  indices 

of  different  phases  of  the  same  soil,  the  semidisperse  model  allows  the 

index  of  refraction  to  be  calculated  for  a given  soil  or  particle  if 

the  dielectric  properties  and  volume  fractions  of  the  constituents  are 

known . 

for  the  KliPii  code,  a more  complete  method  for  specifying  the  index 
of  refraction  for  dust  particles  is  required.  There  are  a number  of 
possible  approaches.  The  simplest  and  most  general  (and  the  easiest 
to  implement  into  the  code)  is  to  allow  the  user  to  specify  as  input 
the  index  of  refraction  for  each  input  frequency.  This  input  option 
is  the  most  accurate  provided,  naturally,  that  the  user  knows  the  index  1 

of  refraction  of  the  dust  particles  (not  the  soil  index)  for  his  case 
of  interest.  In  many  cases  the  user  may  know  various  physical  or  . 

electrical  characteristics  of  the  soils  but  not  the  indices  of  the 
individual  particles.  Later  in  this  section,  we  present  analytic 

methods  of  calculating  the  dust  particle  index  of  refraction  knowing  ] 


either  the  characteristics  of  the  dust  particle  or  of  the  bulk  soil. 

Another  option  is  to  expand  the  built-in  choices  of  soils.  The 
present  choices  of  three  soils  with  widely  spaced  imaginary  indices 
can  be  expanded  to,  say,  ten  soils  with  more  closely  spaced  imaginary 
indices.  The  user  would  choose  one  of  the  ten  typical  soil  types,  and 
the  built-in  index  of  refraction  list  would  consist  of  the  indices 
of  the  dust  particles  formed  from  the  chosen  soil.  Moreover,  the  indices 
would  he  frequency  dependent. 

IVc  now  present  the  semidisperse  model.  Ke  begin  by  defining  the 
complex  permittivity  of  a medium.  The  permittivity  appears  in  the 
literature  in  various  forms  depending  on  how  the  permi tt ivity  was 
measured  or  the  physics  in  which  the  permittivity  appears.  The  most 
common  form  is  the  relative  complex  permittivity: 


= k - i 


k(  1 - i tan  6) 


where 


= relative  complex  permittivity 

= c/e 

0 - 1 
= complex  permittivity  ( 1-  m ) 

= relative  dielectric  constant  (often  reported  simply 
as  the  dielectric  constant) 

= dielectric  loss  factor 

— dielectric  loss  angle  ^ruuy  . oomci. line's  Lne  uiexecci'ic 

phase  angle  0 is  given,  where  0 = tt/ 2 - <5 

= angular  frequency  of  incident  radiation  (rad  s *) 

. • , .....  , . -1. 


material  conductivity  (mho  m ) 
permittivity  of  free  space  (8.85 
/-I . 
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The  complex  index  ol'  refraction  is  related  to  the  relative  complex  per- 
mit t ivitv  bv 


m = v t ' 


Writing  m in  terms  of  its  real  and  imaginary  parts  as 


m = - nn  j 


T (V*  * (ST  ' *) 

y + tan"  6 - l) 


1 + -~T~T~2  " 1 

L K”u)” 
o 


iST  * ;) 


| ( vf  + tan~S  + l) 
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for  a monochromatic  plane  electromagnetic  wave  propagating  in  the  z 
direction  in  a uniform  medium  with  complex  index  of  refraction  m,  the 
electric  field  intensity  can  be  written  as 


:(z,t)  = h eiw(t  - mz/c) 

' o 

-wm,/c  z -iwm-./c  z 
iwt  I R 

= h e e e 

o 

iait  -otz  -i6z 
= Loe  c e , 


where 


) = t + ip  = i - m = propagat  ion  factor  (in  I 

ik,|"l  . . -1 

t = = attenuation  factor  (.nepers  in  I 

e 

3 = ) hi  = phase  factor  (rad  m * j . 

a is  tlie  attenuation  coefficient  for  the  amplitude.  Our 
normally  in  the  power  attenuation  coefficient,  which  is 
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interest  is 


k = da 


- 1 


in 


20w  m 
c in  1 0 


dB  m 


A dust  cloud  consists  of  discrete  dust  particles  dispersed  in  air. 

In  the  general  case,  the  particles  both  absorb  and  scatter  the  incident 
radiation.  In  general,  then,  the  cloud  cannot  be  considered  a uniform 
medium;  and  the  attenuation  coefficient  must  he  calculated  using  the 
Mie  theory  rather  than  the  simple  uniform  medium  expression  of  l.qua- 
tion  72.  Under  certain  restrictive  conditions,  however,  the  cloud 
can  be  considered  a uniform  medium,  and  the  simpler  attenuation  (and 
phase  delay)  equations  can  be  used.  If  the  long  wavelength  limit  ob- 
tains (radiation  wavelength  much  larger  than  the  dust  particle  sites) 
and  absorption  dominates  scattering,  then  the  attenuation  is  proportional 
to  the  dust  cloud  mass  density  and  does  not  depend  upon  the  details  of 
the  particle  sites.  In  tins  simple  case,  a dust  cloud  L:  Jex  of  refrac- 
tion can  be  defined  and  the  uniform  medium  relations  used.  The  exact 
Mie  relations  can  also  be  used  for  this  simple  case,  of  course;  the 
same  attenuation  coefficient  results  from  either  calculation.  When 
scattering  is  significant,  or  the  long  wavelength  limit  does  not  obtain, 
only  the  Mie  relations  can  be  used  for  the  attenuation.  For  backscatter 
calculations,  the  Mie  relations  must  be  used  in  all  cases. 


Consider  a two-phase  mixture,  consisting  of  a disperse  phase  (dis- 
crete elements)  contained  in  a continuous  medium,  as  shown  in  I igurc  3. 
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Figure  3.  A two-phase  mixture. 


Let 


■,*  = relative  complex  permittivity  of  the  two-phase  mixture 
e*  = relative  complex  permittivity  of  the  continuous  medium 
c*  = relative  complex  permittivity  of  the  disperse  elements 
t P = volume  fraction  of  the  disperse  phase 


volume  of  disperse  elements 
total  volume 


Then  the  mixing  rule  for  this  two-phase  system  is  given  by  Hanai , Refer- 
ence (>,  as 


■ * _ t-  * 

'm  “d 


1/5 

(>)  . 1 - 0 


(75) 


Kobschall  (Reference  5J  shows  that  the  permittivity  of  a multi -phase  sys- 
tem can  be  calculated  hy  the  appropriate  repeated  applications  of  this 
two-phase  mixing  rule.  Let 


C = 


e* 

m 


(74) 


I)  - 
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Cubing  both  sides  of  Liquation  73  and  gathering  terms,  the  two-phase 
mixing  rule  can  he  written  is  a cubic  equation: 


. or 
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The  first  form  is  appropriate  if  c*  and  i*  are  known  and  is  the  un- 
known. The  second  form  is  for  L*  and  £*.  known,  with  t*  unknown.  These 

d 111 

cubic  equations  can  be  solved  by  the  standard  techniques,  and  a computer 
routine  has  been  written  to  soLvc  their  tgiven  in  Appendix  C) . 


Our  first  application  of  the  mixing  rule  will  be  the  calculation  of 
the  permittivity  of  a particle  knowing  the  constituent  properties.  The 
particle  will  be  considered  to  be  a mixture  of  solid  material,  water,  and 
air.  first  consider  a large  particle  such  as  a rock,  for  most  rocks, 
conduction  is  electrolytic,  the  conduction  medium  being  an  aqueous  solu- 
tion of  common  salts,  distributed  in  a complicated  manner  through  the 
pore  structure  of  a rock  (Reference  7).  As  long  as  there  is  a continuous 
film  of  water  over  all  the  surfaces  of  the  rock  structure,  the  rock  con- 
ductivity is  a smooth  function  of  the  water  volume  and  water  conductivity. 
The  solid  (essentially  nonconducting)  material  of  the  rock  can  be  con- 
sidered to  be  dispersed  in  a conducting  water  medium.  If  the  water  con- 
tent of  the  rock  is  below  the  cr'tical  saturation  (the  water  film  is  no 
longer  continuous),  then  we  would  have  discrete  conducting  water  elements 
dispersed  in  a nonconducting  medium;  the  conductivity  of  t he  rock  would 
be  greatly  reduced.  Besides  solid  and  water,  a rock  has  air  in  its  pore 
structure.  The  fraction  of  air  in  a rock  is  known  as  its  porosity.  Ke 
will  assume  that  a dust  particle  is  constructed  similarly  to  a rock, 
further,  we  assume  the  water  film  in  the  particle  is  unbroken.  for  the 
| particle,  let 


f = volume  fraction  of  solid  material 
s 


EL* 
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f = volume  fraction  of  air  (porosity) 

f = volume  fraction  of  water 
w 

€* , c* , e*  = complex  relative  permittivities  of  solid,  air,  and 
s u w 


water,  respectively 


Note  that  f + f + f =1 
s a w 


We  calculate  the  permittivity  of  the  particle  by  two  successive  applica- 
tions of  the  mixing  rule,  liquation  76  . First  we  find  the  permittivity 
of  the  solid-water  mixture  by  considering  that  the  solid  is  dispersed 
in  the  continuous  water.  Thus 

= e* 
d s 


‘P  ~ f~  + y 

S Vv 

for  the  first  application.  Let  e*  be  the  solution  for  the  solid-water 
mixture.  The  particle  permittivity  is  found  by  a second  application 
considering  that  air  is  dispersed  in  the  continuous  solid-water  mixture. 
We  use  Lquation  76  again  with 


lor  frequencies  below  about  1 GHz  the  relative  dielectric  constant 
of  water  is  essentially  constant,  and  is  given  by 

K = 87.8  - 0.57T,  (78) 

w 

where  T is  the  temperature  (°C).  It  happens  that  varies  little  with 
soil  type  and  is  approximately  equal  to  5.5  (Reference  5).  We  will 
calculate  the  dielectric  constant  and  conductivity  of  the  particle 
assuming  the  following  parameters: 


39 


79. 1 


Tv 


r 


l 

[ 

i 


»rwpwiP 


K =1.0 
a 


o=() 

a 


5.5 


o=0 

s 


Recall  that 
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where  f is  the  frequency  (Hz).  Applying  the  mixing  rule  twice  as  dis- 
cussed above  leads  to  the  results  shown  in  Table  2.  where  is  the 

-1  P 

conductivity  of  the  particle  (mho  m ). 


Table  2.  Dielectric  constant  and  conductivity  of  a particle. 
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Analytic  fits  to  these  results  are  given  by 

(80) 

(81) 

The  error  in  the  conductivity  fit  is  less  than  1.5  percent.  The  maximum 
error  in  the  dielectric  constant  fit  is  9 percent;  for  all  but  two  table 
entries,  the  error  is  less  than  4 percent.  The  normal  ranges  for  porosity 


Kp  = 50.5  (1.475x2  + l.OlOx  + 0.126)  (0.4055f2  - 1.588^  +1.00) 

o = 0.555O  [2x2  + 1.045x  - 0.025) (0.425f2  - 1 .5f  +1)  mho  m'1 
p w v a a 
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in  rocks  (f  ) is  from  0 to  about  60  percent,  with  a median  value  of  about 

J -] 

15  percent.  For  connate  water  conductivity,  0.05  mho  m is  a low  con- 

ducitvity,  0.3  is  a typical  medium  value,  and  20  is  a high  value  at  radio 

broadcast  frequencies. 

The  DICF  THROW  event  was  an  HL  (high  explosive)  test,  conducted  near 
the  Giant  Patriot  site  on  the  White  Sands  Missile  Range,  6 October  1976 
(Reference  8).  The  charge  for  this  test  was  composed  of  approximately 
628  tons  (570  metric  tons)  of  ammonium  nitrate  fuel  oil  (ANFO) . A com- 
plement of  33  experimenters  and  support  agencies  participated  in  this 
nuclear  simulation  test  event.  SR]  (References  8,9,10)  fielded  a Lilli'/ 
microwave  transmission  experiment  to  measure  the  effects  of  dust  and 
debris  on  signals  passing  through  the  cloud  generated  by  the  blast. 

Botli  amplitude  and  phase  shift  were  measured.  The  measurement  frequen- 
cies are  shown  in  Table  3.  Several  samples  of  loose  crater  material 
were  collected  after  the  detonation  and  analyzed  to  determine  their 
dielectric  properties  and  mass  densities.  Using  the  laboratory  measure- 
ments we  will  compute  the  dielectric  properties  of  the  individual  dust 
grains  by  means  of  the  mixing  rule.  Using  the  dust  grain  values,  we  will 
then  compute  the  propagation  properties  of  the  dust  cloud  by  two  methods: 
first,  by  the  simpler  mixing  rule  method,  and  second  by  the  full  Mie  cal- 
culation. Finally,  we  will  compare  our  calculated  values  of  the  dust 
cloud  attenuation  with  the  experimentally  measured  values. 

We  expect  our  calculated  dust  cloud  attenuations  to  be  equal  or  less 
than  the  measured  attenuations.  This  is  because  the  actual  cloud  has 
not  only  dust  but  the  combustion  products  of  the  628  tons  of  ANFO.  If 
the  extra  attenuation  due  to  these  combustion  products  is  unimportant  at 
a given  frequency,  then  the  calculated  and  measured  attenuations  should 
be  similar.  If  the  combustion  products  attenuation  is  significant,  then 
the  dust  attenuation  will  fall  significantly  below  the  measured  attenua- 
tion. In  the  case  of  a nuclear  device,  the  attenuation  due  to  the  device 
products  will  be  negligible  compared  to  the  dust,  since  the  nuclear  pro- 
ducts represent  an  infinitesimal  fraction  of  the  nuclear-lofted  dust. 
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After  the  test,  soil  samples  were  taken  from  the  resultant  crater. 
There  appeared  to  be  three  characteristic  soils  comprising  the  crater, 
descriptively  labeled  in  Reference  10  as  sand.  Caliche  A,  and  Caliche  B. 
Laboratory  measurements  on  eacli  of  these  soil  types  are  summarized  in 
Table  4 (Reference  10) . 

The  relative  permittivity  of  the  soil  in  the  lab  fixture  is 


- ■ r,  c °oi 

<r-  - i K,-  tan  o + — — = e,.  - 

t f o c u>  t 

L o 


£*  = K 


where  the  subscript  f denotes  values  in  the  lab  fixture.  Let  us  assume 
a linear  variation  of  Kf  with  frequency.  Then  for  the  8 transmission 
frequencies,  the  lab  permittivities  are  given  in  Table  5. 

Next  we  apply  the  mixing  rule  to  the  lab  permittivities  to  calculate 
the  permittivities  of  the  individual  dust  grains.  In  this  case,  the 
mixture  permittivities  are  known  and  we  want  to  solve  for  one  of  the 
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Density  Measurements 


Electrical  Measurements 
of  Soil  in  Lab  Fixture 


Sample 


Density  in  Soil  Grain  Dielectric  Constant  tan  6 
Lab  Fixture  Density  1 ' 


(g  cm'3)  (g  cm'3) 


1.521 


10  GHz 


Cal iche  A 
Caliche  B 


1.207 

1.472 


(mho  m"3 ) 


0.025  0.0024 
0.024  0.0031 


0.021  0.0014 


constituents,  so  liquation  77  is  used.  In  the  lab  fixture  we  consider 
the  soil  mixture  to  be  composed  of  dust  grains  and  air.  The  grains 
are  conducting  and  in  contact  with  each  other,  so  that  the  grains  are 
the  continuous  medium  and  the  air  is  the  disperse  phase.  Thus  in 


iquat i on 


= c*  (grain  permittivity) 


1.0  - i (0)  (air  permittivity) 


e*.  (lab  permittivities) 


air  volume  _ lab  soil  density 
total  volume  grain  density 


0.4059  (sand) 


= <0.5428  (Caliche  A) 


0.4424  (Caliche  B) 


The  results  for  the  dust  grain  permittivities  and  indices  of  refraction 
are  given  in  Table  O. 


Table  5.  Lab  permittivities  for  the  transmission  frequencies. 


Frequency  (MHz) 

378 

608 

413. 

028 

424. 

501 

447. 

447 

1 

273 

503 

2 

891. 

196 

8 

914. 

512 

10 

188. 

024 

Lab  Permittivities 


Caliche  A 


Caliche  B 


2.93  0.187  2.83  0.215 

2.93  0.178  2.83  0.203 

2.93  0.175  2.83  0.199 

2.92  0.169  2.82  0.192 

2.89  0.106  2.79  0.111 

2.82  0.0854  2.72  0.0846 

2.55  0.0686  2.45  0.0651 

2.49  0.0665  2.39  0.0628 


Table  6.  Dust  grain  permittivities  and 
indices  of  refraction. 


M 

' f 

i 

f 

0 

215 

3 

02 

0 

203 

3 

02 

0 

199 

3 

02 

0 

192 

3 

02 

0 

111 

2 

99 

0. 

0846 

2 

94 

0 

0651 

2 

74 

0 

0628 

2 

69 

Frequency 
(MHz  ) 


378.608 

4:3.028 

424.801 

447.447 

1 273.503 

2 o91.196 
8 914.812 

lu  1HK.U24 


Ca 1 iche  A 


0.395  2.20 

0.376  2.20 
0.370  2.20 
0.357  2.20 

0.224  2.18  0.0 

0.180  2.15 

0.143  2.01 

0.139  1 . 98 


Cal  iche  B 
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0.665  2.49  0.133  5.40  0.302  2.32  0 

0.628  2.4Q  0.126  5.40  0.268  2.32  0 

0.616  2.49  0.124  5.40  0.286  2.32  0 

.17  0.594  2.49  0.119  5.40  0.279  2.32  0 

.07  0.343  2.46  0.0695  5.33  0.142  2.31  0 

.86  0.261  2.42  0.0538  5.21  0.163  2.26  0 

.0-3  0.198  2 .24  0.0442  4.75  0.139  2.1; 

.85  0.191  2.20  0.0433  4.63  0.136  2.15  8 


Using  the  dust  gra i n indices  of  Table  • >,  we  next  calculate  the 
propagation  properties  ol'  the  IMfl.  THROW  dust  cloud.  first  we  assume 
that  the  long  wavelength  limit  obtains  and  that  scattering  is  insigni- 
ficant compared  to  absorption.  We  can  then  use  the1  mixing  rule  to  calcu- 
late the  index  of  refraction  of  the  cloud  and  use  !'.|uations  7 1 and  ~J 
to  find  the  phase  shift  and  attenuation  coefficient.  he  use  the  mixing 
rule  of  filiation  To  with  the  dust  grains  dispersed  in  air.  lake  the 
dust  density  in  the  cloud  to  be  10  “ g cm  ' ; this  is  the  estimated  value 
of  the  dust  density  when  the  cloud  had  a diameter  of  about  100  m. 

Table  7 gives  the  permittivities  and  indices  of  refraction  of  the 
dust  cloud.  liquations  71  and  7’J  are  now  used  to  calculate  the  propagation 
properties  of  the  cloud.  We  convert  the  linear  coefficients  into  mass 
penetrated  coefficients  by  dividing  by  the  cloud  dust  density.  These 
mass  penetrated  values  can  then  he  used  for  any  cloud  density.  Table 
8 shows  the  results. 

As  an  example  of  the  use  of  the  propagation  parameters,  we  consider 

a 100  m ray  path  through  a cloud  of  average  dust  density  of  10  “ g cm  . 

(,  . ■> 

lhis  corresponds  to  a mass  penetrated  value  of  10  g m ".  If  the  dust 
were  composed  of  sand,  then  the  dust- induced  phase  shift  and  attenua- 
tion would  he  J.78  rad  and  1.05  dB  at  378. 008  MHz  and  65.)  rad  and 
13.3  dB  at  10188.  (>J4  MHz. 

Next  we  use  the  full  Mie  calculation.  We  use  the  dust  grain  indices 
of  refraction  of  Table  (>  and  the  Wl.l’ll  power- law  size  distribution  de- 
scribed in  the  previous  section.  The  exact  power-law  exponent  for  the 
THROW  dust  cloud  is  not  known,  but  is  thought  to  lie  between  the 
values  1 and  3.  In  order  to  prevent  damage  to  the  experimental  equip- 
ment from  large  soil  missiles  thrown  out  by  the  blast,  an  annular  ring 
of  '>;  soil  around  the  explosive  was  removed  and  replaced  with  sand. 

1 he  dust  cloud  particles  are  biased  toward  small  particles.  The  Mie  ex- 
tinction and  absorption  coefficients  arc  shown  in  Table  9 for  power-law 
exponents  o f 3.5  ("hard  rock"),  4.0  ("unconsolidated  soils"),  and  4.5 
("fine  soils").  The  Mie  scattering  coefficient  is  the  difference  be- 
tween tlie  extinction  and  absorption  coefficients.  At  the  lower 


Table  8.  Dust  cloud  propagation  properties  from  mixing  rule. 


1 


f requeue  tes  and  smaller  particle  distributions  (larger  p),  the  extinc- 
t j on  is  due  entirely  to  the  absorption  and  is  the  same  (within  our  coin 
put.it  ional  accuracy)  as  that  of  table  S calculated  hy  the  simple  mixing 

■i 

rule.  At  the  higher  t requeue ies  and  larger  particle  distributions 
(smaller  p / , the  dust  particles  are  no  longer  small  compared  to  the 
wavelength.  Scattering  becomes  more  important  and  dominates  the  ex- 
tinction. Serious  errors  would  result  it'  the  simple  mixing  rule  were 
used  to  calculate  the  extinction  for  these  eases. 

In  the  HUT.  THROW  dust  cloud  transmission  experiments,  both  phase 

■i 

shift  and  extinction  were  measured,  figure  4 from  Reference  10  shows  1 

the  phase  shifts  on  transmission  path  1 as  a function  of  time  for  the 

three  lowest  frequencies.  We  will  assume  that  the  cloud  is  composed  1 

primarily  of  sand  particles  with  a power-law  size  distribution  exponent 
of  4.5.  Using  the  phase  shift  parameters  from  Table  8 and  the  phase 
shifts  of  figure  -1,  we  can  calculate  the  mean  mass  penetrated  as  a func- 
tion of  time.  We  ignore  the  scintillations  about  the  mean  phase  shift 
1 me . 

figure  f>  shows  the  mean  mass  penetrated.  We  use  the  mass  penetrated 
time  history  and  the  Mie  extinction  coefficients  for  sand  and  p = 4.5  of 
table  d to  calculate  the  mean  extinction  as  a function  of  time  and  fre- 
quency. ligure  (>  shows  the  measured  attenuations  (.Reference  10.1  and  our 
calculated  values.  I he  data  show  strong  scintillations  about  the  mean 
absorption.  Our  calculated  mean  values  show  reasonable  agreement  at  the 
lowest  and  highest  frequencies,  but  do  not  produce  enough  attenuation  at 
the  intermediate  frequencies.  As  we  speculated  earlier,  the  missing  at- 
tenuation at  the  intermediate  freqt  nicies  may  he  due  t’o  the  combustion 
products  of  the  t>Js  tons  of  AMO  in  the  cloud  which  wc  have  neglected. 

Another  possibility  i > that  the  size  distribution  of  the  dust  particles 
may  not  hi  well  represented  hy  a p - 4.5  power  law  distribution.  Be- 
cause part  of  tiie  < soil  was  removed  and  replaced  with  sand,  the 

result  in;  dust  particle  size  distribution  may  he  abnormal. 
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SEE  f I ON  5 

DUST  MODELS  FOR  THE  STEM 
AND  PRECDRSOR  PEDESTAL  REGIONS 


Theoretical  Investigations  of  dust  in  a nuclear  environment  have 
been  carried  out  for  many  years  at  many  different  companies , including 
SA1  (Science  Applications  Inc.),  1ITR1  (Illinois  Institute  of  Technology 
Research  Institute),  ATI  (Applied  Theory  Inc)  and  SJ  (Systems,  Science, 
and  Software)  (see  References  11  to  IS).  The  mechanisms  of  dust  genera- 
tion and  transport  have  been  studied.  Dust  generation  includes  crater 
ejecta,  dust  "popcorned"  from  the  surface  due  to  the  thermal  pulse  of 
the  device  (with  a minor  contribution  due  to  the  neutron  output),  and 
the  subsequent  "scouring"  of  additional  dust  from  the  surface  due  to 
the  outgoing  shockwave  and  the  inrushing  afterwinds.  Transport  involves 
following  the  motion  of  the  dust  after  generation,  both  within  and  out- 
side the  fireball. 

The  general  technique  in  the  transport  studies  is  the  following: 

1.  Assume  an  initial  dust  mass  distribution  in  space  and  a 
particle  size  distribution.  The  initial  mass  distribution  can 
come  from  the  dust  generation  results  or  from  empirical  models. 

2.  Obtain  a model  for  the  air  density  and  velocity  flow  fields 
as  a function  of  time.  Both  large  hydrodynamic  codes  (such 

as  HULL  and  SHELL)  and  idealized  models  (such  as  Hill's  spherical 
vortex)  have  been  used  to  specify  the  fields. 

3.  Follow  the  motion  of  a large  number  of  individual  particles  in 
the  flow  fields.  The  number  and  size  of  the* part  ic 1 es  arc 
chosen  so  as  to  give  good  statistics  for  size  distribution  and 
mass  density.  The  particles  are  assumed  not  to  affect  the 
flow  fields.  The  particle  motion  is  followed  under  the  effects 
of  gravity  and  drag  (from  the  flow). 
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Using  results  from  nuclear  test  experiments,  HP  experiments , and 
the  theoretical  studies,  (lli-TPMl’O  developed  a systems  model  for  the  dust 
within  the  fireball  (Reference  IS).  This  model  was  used  in  the  RAM.' 
and  Will'll  nuclear  weapon  effects  codes,  and  was  subsequently  adopted  for 
the  ROSCOh  systems  code.  Besides  the  fireball  region,  the  stem  and 
pedestal  regions  also  contain  dust.  figure  7 shows  an  artist's  sketch 
of  a surface  or  near-surface  nuclear  detonation.  There  now  exist  sufficient 
data,  both  theoretical  and  experimental,  for  the  development  of  a systems 
model  for  the  dust  within  the  stem  and  pedestal  regions.  In  this  section 
we  develop  such  a model. 


Figure  7.  Sketch  of  the  nuclear  dust  regions. 


I 
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We  first  consider  the  dust  pedestal  region.  A typical  pedestal 
region  will  extend  out  to  about  three  initial  fireball  radii  and  be 
SO  to  100  meters  thick,  containing  roughly  0.03  Ml'  of  dust  per  MT  of 
yield  (Reference  1 1 J . The  formation  of  a large  dust  pedestal  region 
is  a unique  characteristic  of  a nuclear  detonation.  About  40  percent 
of  the  total  yield  of  a nuclear  detonation  is  emitted  from  the  fire- 
ball in  the  form  of  thermal  radiation.  The  soil  surface  around  the 
fireball  region  is  irradiated  by  this  thermal  flux.  The  weapon  neutron 
and  X-ray  output  may  also  irradiate  the  soil  surface;  but  the  range  of 
these  outputs  is  so  short  that  the  affected  soil  is  within  the  fireball 
or  stem  regions,  and  so  they  do  not  contribute  to  pedestal  function. 

The  thermal  irradiation  causes  the  top  0.3  to  1 mm  layer  of  soil  to  be 
thrown  off,  forming  a low-lying  dust  layer.  For  wet  soils,  steam  pro- 
duction is  probably  the  throw-off  mechanism.  For  dry  soils,  a "pop- 
corning"  effect  is  observed,  probably  due  to  the  release  of  hydrated 
water  (Reference  12J. 


The  air  in  this  low-lying  dust  region  is  heated  by  a variety  of 
mechanisms.  There  is  conduction  from  the  hot  soil  surface,  convection 
within  the  layer,  steam  production  in  the  case  of  wet  soils,  and  con- 
duction heat  from  the  hot  thrown-up  soil  particles.  Once  the  dust 
layer  is  formed,  the  dust  particles  absorb  the  incident  thermal  radia- 
tion and  efficiently  transfer  the  absorbed  heat  to  the  air  by  conduc- 
tion. The  net  result  is  the  formation  of  a thin,  hot,  dusty  air  layer 
above  the  surface.  Sound  propagates  faster  in  heated  air.  As  the  blast- 
produced  shock  wave  propagates  outward  from  the  burst,  that  portion  of 
the  shock  propagating  in  the  heated  region  travels  faster  than  the  main 
shock;  it  becomes  the  shock  precursor.  Figure  8 from  Reference  12  shows 
a sketch  of  the  shock  behavior.  The  outrunning  shock  wave  scours  some 
additional  dust  from  the  soil  surface.  The  formation  of  a shock  pre- 
cursor causes  large  vertical  air  velocities  behind  the  shock  front  which 
carry  the  dust  upward  50  to  100  meters. 


The  dust  production  from  a 
cantly  from  that  of  a nuclear 
hot  air  layer  existing  before 


conventional  explosion  differs  signifi- 
explosion.  Conventionally,  there  is  no 
the  shock  arrival.  The  outrunning  shock 
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rsor  development. 
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can  scour  dust  from  the  surface,  hut  there  are  no  precursor- induced 
vertical  velocity  components  to  carry  the  dust  to  relatively  high  alti- 
tudes. Thus  the  total  dust  production  is  proportionally  less  and  the 
dust  is  confined  to  the  near  surface  layer.  The  nuclear-producod  dust 
is  truly  a "precursor"  pedestal. 

Tor  the  development  of  our  systems  model  for  the  dust  pedestal, 
we  rely  primarily  upon  the  work  reported  in  References  11,  12,  and  13. 

In  these  reports  the  authors  have  studied  the  pedestal  formation  and 
time  history  from  both  a theoretical  and  experimental  approach.  The 
experimental  data  consists  primarily  of  photographic  records  of  nuclear 
test  events  at  the  Nevada  Test  Site;  there  arc  also  a few  dust  density 
measurement s . 

We  begin  with  the  geometrical  modeling  of  the  dust  pedestal.  from 
test  observations,  the  thermally  produced  dust  layer  is  rapidly  lofted 
to  its  maximum  altitude  after  passage  of  the  shock  wave.  The  radius 
grows  out  to  a fairly  well  defined  maximum  radius  (where  the  shock  be- 
comes too  weak  to  loft  the  dust),  figure  9 from  Reference  12  shows 
the  time  history  of  some  typical  dust  profiles  predicted  by  the  I’kklUJM 
code.  SAI  approximates  the  final  pedestal  dust  cloud  by  the  simple 
geometry  shown  in  figure  10.  Utilizing  the  photographic  data,  an 
empirical  model  for  R^  and  was  developed  in  Reference  12.  This 

model  is  shown  in  figure  11;  also  shown  is  an  estimate  of  IT,,,,  developed 
from  I’RfDUM  calculations.  I'or  our  pedestal  model  we  choose  the  simple 
geometry  of  an  annular  cylinder;  the  inner  radius  is  taken  as  the  stem 
radius.  Our  geometry  is  shown  in  figure  12. 

As  the  burst  altitude  is  increased  above  the  ground  surface,  the 
thermal  flux  at  the  ground  and  the  strength  of  the  outrunning  shock 
wave  become  weaker.  Tventually  an  altitude  is  reached  at  which  no  pre- 
cursor is  formed.  Reference  12  estimates  that  this  point  is  readied  at 
an  altitude  of  about 

llNlu(ftJ  = 075  W1A>  (KT)  . (kku 
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Figure  9.  Shot  PRISCILLA— sweep-up  dust  laye1 
at  various  times. 
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Figure  10.  PREDUM  [predictions  for  shot  DOG  and  simplified 
geometrical  fit. 

The  units  customarily  used  in  the  stem  and  [Pedestal  references  are  feet 
for  distance  and  kilotons  for  yield.  I he  working  units  of  most  weapons 
effects  codes  such  as  Mi’ll  are  kilometers  for  distance  and  megatons  for 
yield.  We  will  general ly  write  the  same  equation  twice,  once  in  each 
set  of  units.  The  distance  and  yield  units  wall  he  explicitly  indicated 
within  the  equation;  the  units  of  the  other  variables  which  do  not  change 
will  not  he  explicitly  indicated  in  the  equation,  hut  will  he  specified 
elsewhere.  Thus  the  maximum  altitude  above  which  no  precursor  is 
formed  is  also  written  as 


1 MAX 


I 

(4  .S 


(1.71  x 10  ♦ 1.00  A,1/5(MT) 


0 < nfi(kni)  <_  0.62SStl/3(Mi; 


US  X 10'5Ub(UO  * 2.S7  W1/3(MT)  0.fc25W1/,i  [MT)  < Hgfkin)  <_  2 . 06W1/3 (MT) 

(H4hj 

figure  13  from  Reference  12  shows  the  estimates  for  the  height  of 
the  final  pedestal  dust  cloud;  shown  arc  the  test  data  points,  the 
"host  estimate"  (empirical  model),  and  the  I’RfDUM  code  predictions. 

IV e adopt  the  following  fit  (also  shown  on  figure  15): 
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he  will  model  the  time  history  of  the  pedestal  radius  as  it  grows 
from  the  stem  radius  to  its  max i muni  radius.  he  ignore  the  short  build- 
up time  that  it  takes  t lie  dust  to  reach  its  maximum  altitude,  and  take 
the  pedestal  height  as  a constant  equal  to  lip  Consider  a burst  at 

an  altitude  low  enough  for  a precursor  to  he  formed.  The  shock  wave 
travels  downward  from  the  burst  point,  intersects  the  ground,  and  then 
travels  outward  along  the  ground  surface.  The  precursor  forms  and  the 
dust  pedestal  is  generated,  he  take  as  the  pedestal  outer  radius  time 
h i story , 
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t^u  = time  of  arrival  of  the  shock  at  the  ground  (s) 
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Figure  13.  Dust  height  predictions  (PREDUM)  compared  to 
empirical  dust  model  and  data. 

1 = ground  range  of  the  shock  wave  (ft  and  km  J 

t . = time  at  which  the  precursor  forms  (sj 

t = time  constant  which  determines  the  rate  at  which 

e 

the  final  radius  is  reached  (s). 

In  our  model  there  is  no  dust  pedestal  until 


Rp(tl  ’ «slt) 


where  R^.(t)  is  the  radius  of  the  stem  and  the  inner  radius  of  the  pedestal 
(.R^ltj  will  he  given  latcrj  . The  equation  for  R (tl  for  t _ t_  v is  of 
the  form  given  in  Reference  13. 
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The  slant  range  from  the  burst  point  to  the  shock  front  i' 
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l'efinc  a scaled  slant  ranee  as 
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I hen  the  time  in  seconds  at  which  the  shock  front  reaches  ground  radio- 

1<  . , , is  riven  in  Reference  la  as 
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where 


f = lo'V'  -'i  KI  i - Ilf  IV  ' I Ml  ; 


This  expression  is  a fit  to  tlie  shock  rau'ius  history  -is  calculated  hy  a 
shock  hydrocode.  I he  t i nv  of  arrival  of  tile  shock  at  the  ground  is 
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he  can  invert  liquation  89  to  find  the  shock  front  rail  i us  as  a i 
of  time.  Define  a scaled  time 


Then  inverting  liquation  89, 
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Reference  15  gives  the  time  constant  t as 
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I igure  14  from  Reference  12  shows  the  average  dust  density  in  the 
final  dust  cloud  as  a function  of  yield  and  altitude,  as  calculated  hy 
the  PRhOUM  code.  An  excellent  fit  to  the  data  is 
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Figure  14.  Average  dust  density  in  the  final 
dust  pedestal . 

What  is  the  size  distribution  of  the  particles  in  the  dust  pedestal? 
We  arc  assuming  that  the  particles  came  primarily  from  the  to]'  0.3  to  i 
mm  of  t he  soil  surface,  so  the  particles  arc  small.  We  will  assume  that 
the  pedestal  particles  have  the  same  size  distribution  as  the  fireball 
dust  particles,  but  that  the  maximum  size  particle  in  the  pedestal  is 
1 nun  (,0.1  cm).  For  evaluation  purposes  the  fireball  dust  size  distribu- 
tion is  divided  into  eight  size  groups  spanning  particle  diameters  from 
0.001  to  10  centimeters.  We  take  the  same  small-size  limit  of  0.001 
centimeters.  For  evaluations,  we  divide  the  pedestal  dust  sizes  into 
the  four  size  groups  given  in  Table  10.  Since  we  are  assuming  the  same 
size  distribution  for  the  pedestal  dust  as  for  the  fireball  dust,  we  can 
use  the  fireball  dust  attenuation  and  backscatter  formulas  for  the 
pedestal  dust  propagation. 


Table  10.  Pedestal  dust  sizes. 


I 


Group 

Minimum  Diameter 
of  Particle  in  Group  (cm) 

Maximum  Diameter 
of  Particle  in  Group  (cm) 

1 

0.0"1 

0.004 

2 

0.004 

0.01 

3 

0.01 

0.04 

4 

0.04 

0.1 

As  In  the  fireball  case,  we  assume  that  each  size  group  is  uniformly 

dstributed  within  its  own  (annularj  eylindrical  volume.  All  'cylinders 

will  be  assumed  to  have  the  same  radial  dimensions  at  all  times,  inner 

radius  R (t  ) and  outer  radius  R (t).  All  cylinders  will  initially  have 
s P 

the  same  altitude,  II,,  ; after  a suitable  delav  interval  the  height  of 

l ’max 

each  cylinder  will  he  allowed  to  fall  with  its  characteristic  velocity. 
We  do  not  have  any  data  for  the  appropriate  delay  time  before  the  dust 
particles  begin  to  settle.  We  arbitrarily  take  the  delay  time  to  be  5 
times  the  e-folding  time  of  U (t)  in  liquation  S6, 


t ..  = t 

ip  pl- 


at 
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ihe  altitude  history  of  each  group  cylinder  is 
‘i’ma.x 


(II  J 
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t < tt. 
- fp 


i 
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Pmax 


v . (t  - t ..) 

pi  1 
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where  v ^ is  the  characteristic  fall  velocity  of  size  group  i. 

The  terminal  velocity  equation  of  a spherical  particle  falling  under 
gravity  in  still  air  is 


Aiih  v- 
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where 

V = terminal  fall  velocity 
p = ambient  air  density 
p = particle  density 
d = diameter  of  spherical  particle 
g = acceleration  of  gravity- 

{•or  a spherical  dust  particle  we  can  use  (Reference  11  j 

CU  - ir  * °'44  • 

e 

where  R is  the  Reynolds  number,  given  by 

p Vd 
R = JL_ 

6 "a 

tiiid  M is  the  air  viscosity.  Thus  the  velocity  equation  becomes 
a 

1 pa  ^ 18Ua 

i rr r * 7 v - B = 0 ■ 

V1 

which  lias  the  solution 


21  3gpd  27y 


This  terminal  velocity  equation  is  more  familiar  in  the  limits  of  large 
and  small  dust  particles.  For  small  dust  particles  (d  < 0.01  cm) 

24  , 

S = IT  and 

e 

v = , (106a) 

18pa 

which  is  Stokes  law.  l or  large  dust  particles  (d  ;>  0.1  cin)  , 


V = \/^V 

V ->,jalT) 


(106b.) 


68 


which  is  the  dynamic  drag  equation.  We  assume  the  following  numerical 
values  : 


p = 1.225  > 10  J g cm  J (sea-level  air) 

p = 1.6  c cm  J 
P _? 

g = 980.6  cm  s 
-4 

p = 1.81  x 10  poise 

With  these  values  the  terminal  velocity  equation  becomes 


^3.9894^  r . 6 . 

: — + 6 . 24. ->8  ■ 10  d 

l d / 


3.9894 


cm  s 


( I<J7  ) 


where  d is  in  cm.  We  assume  that  each  size  group  falls  with  the  terminal 
velocity  characteristic  of  the  smallest  particle  in  the  group.  The  dust 
particles  are  not  spherical  nor  are  they  falling  in  calm  air.  Hence  we 
expect  the  effective  fall  velocity  to  be  less  than  our  simple  calculation. 
Taking  the  smallest  particle  in  the  group  assures  a slower  rate.  'I able 
11  gives  the  fall  velocities  for  each  size  group. 


Table  11.  Fall  velocities  for  pedestal  dust  size  groups. 


Group 

Diameter  of  Smallest 
Particle  in  Group  (cm) 

Fall  Velocity 

cm  s'1 

km  s"1 

ft  s'1 

1 

0.001 

0.782 

k£> 

1 

CNJ 

CO 

r-'v 

2 . 57  ^ 

2 

0.004 

12.4 

1.24-4 

0.408 

3 

0.01 

71.8 

7.18'4 

2.36 

4 

0.04 

410 

4.10'3 

13.4 

We  next  develop  a systems  model  for  the  dust  in  the  nuclear  stem  le- 
gion. The  stem  region  is  more  complicated  to  model  than  the  pedestal  re- 
gion. lor  very- low-a 1 t i tude  bursts,  the  fireball  itself  intersects  the 
ground  and  dust  is  directly  injected  into  the  fireball.  The  fireball  and 
stem  are  connected  and  rise  together;  additional  dust  is  injected  into  the 


fireball  via  the  stem.  As  the  burst  altitude  is  increased,  the  fire- 
ball does  not  reach  the  ground  surface.  A dust  stem  is  formed  .uid 
rises  to  intersect  the  rising  fireball;  dust  is  injected  into  the  tire- 
hall  only  via  the  stem.  At  still  higher  hurst  altitudes,  the  stem  does 
not  reach  the  rising  fireball;  the  fireball  region  is  dust  free.  1 vcn- 
tually  a high  enough  hurst  altitude  is  reached  so  that  no  stem  is  formed. 

We  take  the  geometry  of  the  stem  region  to  he  cy  1 indr  i>--t  1 ; sec  ! i..ure 
12.  Wc  begin  by  modeling  the  stem  radius  time  histor>.  1 igurc  la  fron, 
Reference  11  shows  the  scaled  stem  radius  as  a function  of  scaled  time 
for  a number  of  very  - low -a  It  it  tide  bursts.  These  data  were  taken  from 
photographic  data  of  the  test  events.  We  adopt  the  following  fit  to 
the  best -est imate  curve  shown  on  ! igure  la: 
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There  is 

no  discernible  variation 

of  k,.(tj  with  height  of  burst  for  the 

limited  altitude  range  of  Figure  15.  We  assume  therefore  that  there  is 
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no  functional  dependence  upon  height  of  burst  for  the  whole  narrow 
altitude  range  for  which  a stem  is  formed.  Ke  restrict  K ( t ) to  be  the 
minimum  of  that  given  in  liquation  loti  and  the  shock  front  radius  of 
kquation  '.>(>.  Note  that  for  a burst  altitude  less  than  6.S.9K*''1  feet 

kilometers],  the  stem  radius  eventually  expands  enough  to  en- 
gulf the  entire  dust  pedestal. 

define  a scaled  height  of  burst  of 


where 

Kj  = fireball  pressure  equilibrium  radius  (the  radius  of  the 

fireball  when  the  fireball  internal  pressure  falls  to  the 
ambient  atmospheric  pressure). 

lhon  from  observations  of  nuclear  test  events,  the  following  general  state- 
ments can  be  made  about  the  rise  behavior  of  the  stem  (Reference  1 B)  . lor 
scaled  heights  between  about  2.5  to  5,  there  is  no  mixing  of  the  fireball 
proper  and  the  surface  dust  layer.  The  stem  does  not  begin  rising  until 
~'0  to  50  seconds  after  detonation,  giving  the  larger  particles  initially 
lotted  ample  time  to  fall  back.  The  rising  dust  column  (the  stem)  docs 
not  reach  the  rising  fireball  until  a few  minutes  after  burst,  lor  scaled 
burst  heights  greater  than  about  5,  the  column  of  rising  dust  does  not 
reach  the  rising  fireball  cloud  before  the  cloud  stabilizes  (at  about 
7 minutes) . 

Based  on  a very  limited  set  of  test  data,  we  adopt  the  following  some- 
what arbitrary  model  for  the  time  history  of  the  stem  altitude.  l.et 

11^  (t)  = altitude  of  the  top  of  the  stem  (ft  and  km) 

1 1 But (t ) = altitude  of  the  bottom  of  the  fireball  (ft  and  km) 

t..  _ = time  at  which  the  fireball  stabilizes  (s) 

r do 

t^  = time  at  which  the  rising  stem  intersects  the  fireball 

bottom  (s) 

t = delay  time  before  the  stem  starts  to  rise  (s) . 


12 


! 


L 


Hu-  stem  altitude  will  always  be  limited  by  the  fireball  bottom  altitude; 


“s'* 


lf,,.,(t)  at  all  times 
hU  I 


We  assume  there  is  no  stem  formation  until  the  downward-traveling  shock 
wave  from  the  burst  intersects  the  ground  surface.  We  then  let  the  alti- 
tude rise  with  the  same  velocity  as  t he  outrunning  shock  until  the  stem 
altitudes  reaches  either  the  dust -pedesta 1 maximum  altitude  or  the  fireball 
bottom,  depending  upon  the  burst  altitude.  for  bursts  with  hc„  < 2 , take 
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ksn<-t  1 -•  "box'*  J 

1110) 

'wtJ 

kSH(t)  > l.B0T(t, 

where  is  given  by  fquation  91,  K by  liquation 

9 b , 

and  llB0T  is 

given  hy  the  fireball  phenomenology  model.  Thus  for 

very 

low  bursts, 

we  let  the  stem  rise  until  it 

intersects  the  fireball 

; it 

remains  at  the 

t’i  reball  bottom  from  then  on. 

lor  bursts  with  h^B  2,  we  take 

(° 

1 1 <511 

H.ltj  = Ksl|(tj 

Rsn(tj  - "''max 

(ini  : 

'"'•MAX 

< < <1. 

wliere  llpMAX  is  the  maximum  altitude  reached  by  the  dust  pedestal  and  is 

given  by  fquation  S5 . for 

bursts  above  a scaled  height 

of  2,  we  let 

the  stem  rise  to  the  pedestal 

altitude  and  then  hold 

it  at  that  altitude 

tor  the  delay  time  tj  . Take 

the  delay  time  to  be 

1" 

hSB  1 2 

tn  = r’0  (hSB  2j 

2 i hSB  1 5 • 

(1121 

(so 

hSB  " 3 

A 


Wo  simply  Lot  tho  delay  time  increase  from  0 to  50  seconds  as  the 
scaled  height  of  burst  increases  from  2 to  5. 

After  the  del;  v time  is  over,  the  stem  rises  again.  We  assume  that 
the  stems  for  hursts  with  scaled  heights  between  2 and  5 rise  and  inter- 
sect the  fireball  bottom.  We  take  the  time  of  intersection  as 


1 1 BS  ^ hSB 


lor  bursts  with  sealed  altitudes  between  5 and  0 , we  assume  the  stem  does 
not  reach  the  fireball,  but  stabilizes  at  a lower  altitude  given  by 


lb  - hSBj 


‘il<)Tlt|  BS  ' 


‘‘>  - hSB  - b 


lor  bursts  with  h,.„  • 0 we  assume  no  stem  forms.  As  a comparison,  the 
S 15 

maximum  scaled  burst  height  for  which  a dust  pedestal  is  formed  is 
hs(.  a 2.(i.  The  Wl.l’H  fireball  dust  model  assumes  that  dust  is  lofted  into 
tiie  fireball  for  hursts  with  h ...  5.  Our  stem  model  is  consistent  with 

the  current  fireball  model.  Table  12  gives  the  present  model  dust  region; 
as  a function  of  scaled  burst  height. 


lor  h^B  ' 2,  define  the  average  rise  velocity  \'s  of  the  stem  as  it 
rises  from  the  pedestal  altitude  to  either  the  intersection  with  the 
fireball  bottom  (2  hs(  ^ 5)  or  to  its  maximum  altitude  (5  h^  ;;  (>.i 


llB0T(tlS)  ~ "''MAX 
\ 1 IS  " t|) 


2 - hSB  - 5 


r^MAX  ~ “'MAX 
1 1 I BS  ' t|) 


3 hSB  1 (1 


We  choose  the  simplest  time  rise  history,  namely  a constant  velocity  rise 
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Table  12.  Dust  regions. 


Scaled  Bui  it  Height  (h^) 

Dust  Regions 

0 - 2.6 

Stem,  fireball,  pedestal 

2.6  - 3 

Stem,  fireball 

3 - 6 

Stem 

> 6 

No  dust  regions 

| "'’max  + Vr  ' V 


li>  ::  ! r is 


1(, 


t t 


IS 


sb 


<> 


("‘’max  + Vs  lt  " tnl 
1 is  1 1 i = ' 

I1  'SM.VX 

he  have  expressed  the  rise  model  is  a function  of  sealed  height  of 
hurst.  lor  comparisons  with  the  unsealed  units,  we  can  use  the  conversions 
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Next  we  model  the  dust  densities  in  the  stem  region.  he  consider 
only  that  dust  which  is  lofted  by  the  air  velocity  flow  fields.  We 
ignore  the  crater  ejecta  and  fallback,  lor  scaled  burst  heights  above 
about  0.03,  the  crater  becomes  a compaction  crater  with  very  little 
ejecta,  so  the  lofted  dust  is  entirely  composed  of  material  "popcorned" 
or  scoured  off  the  ground  surface.  lor  bursts  low  enough  to  form  an 
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iTiitiT,  thf  fireball  intersects  the  ground ; any  systems  dust  ef- 
fects would  he  comp  1 et  el  y ot)seured  by  the  fireball  itself.  On  1 y after 
t fie  fireball  has  r;  sen  sufficiently  and  before  the  ejecta  falls  hack 
would  the  stem  electa  be  significant.  I lie  fireball  dust  model  accounts 
tor  the  crater  ejecta  which  is  entrained  into  the  rising  fireball. 

Ue  lose  little  systems  significance  by  ignoring  the  stem  ejecta. 

In  Reference  11  the  early-time  lofted  stem  dust  densities  hare  iio.ii 
est  f mated  In  combining  VORIUIM  v ode  predictions  with  the  available  photo- 
g rapine  data  of  the  stem  time  history.  There  are  fairly  large  uncertainties 
in  the  calculation  of  the  stem  densities.  As  would  be  expected,  there  arc 
no  direct  experimental  measurements  of  the  dust  density  in  a nuclear  s t 
with  which  to  compare  the  calculated  estimates. 

I i gu res  1<>  through  2t>  from  Reference  I I show  the  estimated  dciisit  ies 
for  a range  of  burst  yields  and  altitudes.  The  solid  lines  in  the  fig- 
ures are  the  best  estimates  of  the  densities.  '1  here  is  no  definite  func- 
tional dependence  of  dust  densities  upon  height  of  burst  evident  from 
the  calculations;  therefore  we  assume  no  height  of  burst  dependence.  he 
adopt  the  best  estimate  curves  as  our  primary  data,  and  devise  the 
following  fit.  The  curve  fit  is  also  shown  in  the  figures  for  comparison. 


Let 

/ 5.2  x 10~'V'a’  = 5.2  x Id"- 1\ V'r  K ■■  5 Ml  = 50(10  KT 

I (kl)  (.Mlj 

I (5 . !>  x Iff'  g cm'"'  U n 5 Ml  = 5000  KT 


1 / 

i 1 1 0 1 

Rlltti  = 2-.0lVkT| 

(120a 

R.ikm,1  = 0.-521^^ 

( 1 20b 
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Fiyure  20.  Early  time  stem  densities  for  0.50-MT  burst  at 
SO  W^( KT)-feet  altitude. 
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Figure  25.  Early  time  stem  densities  for  5.0-MT  burst  at  50  W1/,3(KT) 
feet  altitude. 
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Figure  26.  Early  time  stem  densities  for  5.0-MT  burst  at  150  W^(KT) 
feet  altitude. 
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Figure  27.  Early  time  stem  densities  for  20.0-MT  burst  at  50  W^(KT) 
feet  altitude. 
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the  first  term  accounts  for  the  early  time  rapid  falloff  while  the  second 
term  accounts  for  the  late  time  behavior.  1 he-  late  time  functional  de- 
pendence is  equivalent  to  assuming  that  the  total  dust  mass  per  unit  lengtl 
of  the  stem  remains  constant  as  the  stem  radius  expands;  ie,  this  is  a 
oust  mass  cense r\  at  ion  assvunpt  i on  . 

Code  calculations  indicate  that  the  air  flow  fields  in  the  stem 
region  are  strong  enough,  to  loft  soil  particles  of  about  1 centimeter 
or  so  in  diameter,  if  these  particles  are  initially  lofted  into  the  flow 
fields  from  the  soil  surface  (Reference  l-1  ).  l\e  extend  the  ['article 
size  groups  chosen  for  the  pedestal  model  up  to  1 centimeter  in  diameter 
for  the  stem  size  groups;  Table  la  shows  tlu*  stem  size  groups. 


Table  13.  Stem  dust  size  groups. 


Group 

Minimum  Diameter  of 
Particle  in  Group  (cm) 

Maximum  Diameter  of 
Particle  in  Group  (cm) 

1 

0.001 

0.004 

2 

0.004 

0.01 

3 

0.01 

0.04 

4 

0.04 

0.1 

5 

0.1 

0.4 

6 

0.4 

1.0 

In  Reference  la  the  trajectories  of  dust  particles  up  to  1 centi- 
meter diameter  are  studied  by  means  of'  the  SUM  code.  I he  calculated 
trajectories  depend  in  a complex  manner  upon  yield,  height  of  burst,  and 
initial  assumed  particle  injection  velocity.  The  calculations  considered 
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times  out  to  uO  seconds  and  burst  heights  between  J()K^'’(KI|  and  JOOU'* / 
(kll  feet.  These  theoretical  data  are  not  complete  enough  to  allow  an 
analytic  model  oh  stem  particle  motion  to  be  developed;  moreover  we  do 
not  know  the  proper  injection  velocity  to  use.  Hence,  we  will  develop 
a very  simplified  stem  particle  motion  model  for  use  until  more  complete 
theoretical  data  are  available. 


he  begin  by  calculating  the  maximum  altitude,  (11  j that  each  sice 
group  can  be  carried.  IV e take  the  characteristic  air  veioettv  to  be 
that  of  the  average  stem  velocity,  V , of  Iquation  115,  IVc  then  substi- 
tute \ into  the  terminal  velocity  1 quation  1 0 1 and  solve  for  , . ibis 

is  the  densitv  a flow  field  of  veloeite  V must  have  in  order  to  iust 
* ' s 

offset  the  gravity  force  of  a particle  of  diameter  d 


54i’a 

'd.\ 

l : 


cm 


i - . 


he  assume  that  the  stem  flow  field  carried  a particle  of  diameter  u. 

i 

upward  to  an  altitude  (.11^.).  where  the  ambient  density  is  (.  The 

flow  field  can  loft  only  those  particles  for  which 


where  f is  the  ground  level  air  densitv.  lor  the  character i st ic  particle 
diameter  dj  of  Iquation  lid,  take  the  maximum  diameter  of  the  particles  in 
group  i . 

The  altitude  of  each  site  group  cylinder  is  identical  to  the  stem 
altitude,  II..  ( tj,  until  the  altitude  fli^J  . is  reached.  The  group  cylinder 
altitude  stops,  then  falls  as  the  stem  flow  fields  weaken.  he  again  use 
the  minimum  diameter  of  the  group  site  to  calculate  the  fall  velocity 
(.again  to  crudely  compensate  for  the  fact  that  the  air  is  certainly  not 
quiet;.  IVe  use  the  terminal  velcity  equation  to  calculate  two  fall 
velocities  Y^.  and  V^.  , i here  V^.  is  the  velocity  at  altitude  (11^). 
and  V . is  the  velocity  at  ground  level.  for  small  particles  these 


’ll 


velocities  will  lie  about  the  same  and  we  can  use  the  average  as  a '.distant 
fall  rate.  lor  larger  particles  the  fall  rate  will  be  larger  at  the  high 
altitude  and  will  diminish  as  the  particle  encounters  denser  air  at  the 
lower  altitudes.  Assume  that  the  velocity  is  a simple  linear  function  of 
altitude.  then 


v(HJ  - v + rr — (V  - V ) 

o II...  h o 
SM 


Solving  the  first-order  differential  equation  for  II,  we  have 


r V -V 

v, ii . - — °-  it  - t,j 

lilt)  = 11  - — 1 1 - e SM  fs 

1 • SM  V.-V 
h o 


t > tr  U-5 
f s 


where  t^-  is  the  time  the  fall  begins,  lor  small  particles 

V,  ~ Y and  we  take 

h o 


iiuj  = hsm  - vo (t  - tfsj 


r ’ tfs 


lienee  the  t inn.  history  of  the  altitude  of  each  group  cylinder  i: 


l II  It)  t < (t.J. 

nat)  = s b 1 

I lift)  (cqs  Ids  and  I d(i  j t > (t^J^ 


and  ( t ..  ].  is  the  time  at  which  II  ft,.  ) = 111.,,)., 
fs  i s ts  yM  i 
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MIE  COMPUTER  ROUTINES 
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In  this  appendix  we  give  computer  listings  for  two  versions  of  the 
improved  Mie  routine.  The  first  version,  MISCAT,  calculates  the  Mie 
efficiencies  for  extinction,  scatter,  and  backscatter.  The  second 
version.  Mil:,  calculates  the  Mie  efficiencies  for  extinction  and 
scatter  and  the  unnormalized  scattering  pattern.  From  either  version, 
the  Mie  efficiency  for  absorption  can  be  found  from 


^ABS  = ^HXT  ‘ ^SCA 


(A-l) 


The  unnormalized  scattering  pattern  calculated  in  Mil:  is 


V,j 


S.C0J 


(A-2) 


where  S(8)  is  the  scattering  pattern  given  by  liquation  5 in  Section  2. 
For  the  Mil:  version,  the  backscatter  efficiency,  if  desired,  can  be 
calculated  from 


Q 


BKS 


4SyCH) 


2 

oT 


(A- 3) 


where  a is  the  dimensionless  size  parameter  defined  in  Section  2. 

The  complex  function  routine  ANF  evaluates  the  complex  function 
A lY)  (Equation  19  of  Section  2)  and  is  used  by  both  the  MIF.  and  MISCAT 
routines.  Computer  listings  for  ANF  are  also  included  here. 
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SUBROUTINE  M 1 SCAT ( XN,  DR,  01,  DBS*  OEXT.  OSCA  ) 


THIS  IS  A MODIFIED  WOE  ROUTINE  - MODEL  EQUATIONS  DOCUMENTED  IN 

REPORT  GE77TMP-22 

THIS  ROUTINE  USES  MIE  THEORY  TO  CALCULATE  THE  EFFICIENCIES  FOR 
SCATTER (TOTAL) » BACKSCATTER  AND  EXTINCTION  FOP  A SINGLE  UNIFORM 
SPHERICAL  PARTICLE 

INPUTS 

XN  = NORMALIZED  SIZE  PARAMETER,  WHICH  EQUALS  TWO  PI  TIMES 

THE  RADIUS  OF  THE  SPHERE  DIVIDED  BY  THE  WAVELENGTH  OF  THE 
INCIDENT  RADIATION 

DR  = REAL  PART  OF  THE  COMPLEX  INDEX  OF  REFRACTION  OF  THE 

SPHERE 

DI  = IMAGINARY  PART  OF  THF  COMPLEX  INDEX  OF  REFRACTION  OF  THE 

SPHERE 

MOTE  THAT  THF  OPLFX  IMQFX  OF  REFRACTION  IS  ASSUMED  TO  BE 
N = DR  - 1*0 1 


OUTPUTS 


SCATTER(TOTAL)  EFFICIENCY,  WHICH  EQUALS  THE  (TOTAL) 
SCATTERING  CROSS  SECTION  DIVIDED  BY  THE  CROSS  SECTIONAL 
AREA  OF  THE  SPHER£<  SIGMA/ (PI*RADIUS**2) 

BACKSCATTER  EFFICIENCY,  WHICH  EQUALS  THE  SCATTERING  CROSS 
SECTION!  IN  THE  BACKWARDS  DIRECTION)  DIVIDED  BY  THE  CROSS 
SECTIONAL  AREA  OF  THE  SPHEREt  SIGMA/ <PI*RADIUS**2> 
EXTINCTION  EFFICIENCY,  WHICH  EQUALS  THE  TOTAL  (SCATTERING 
♦ ABSORPTION  ) CROSS  SECTION  OF  THE  SPHERE  DIVIDED  BY  THE 
CROSS  SECTIONAL  AREA  of  the  sphere 


DIMENSION  ANR ( 200 ) 

COM°LFX  D , Z » FM 1 , EM2  » EM , ANF , ANR , ANz • AN  »BN ,QBSC 
D = CMPLX ( DR,  -DI  ) 


XsAMlNK  100.,  XN  ) 

Z = X * D 

EMI  = CMPLX ( SIN ( X 
EM2  = CMPLX ( COS(  X 
OBSC  = ( 0,,  0,  ) 
OEXT  = 0, 

OSCAsO.O 
NX  = 1.5  * X 
NXsMAXOf  2,  NX  ) 

ANR (NX)  = ANF<  NX,  i 
NXMI  = NX  - I 
DO  10  I = 1,  NXMI 
N r NX  ♦ 1 - I 
CN  = FLOAT ( N ) 
ANR(N-l)  = Cfl  / Z - 
CONTINUE 
X2  = X **  2 


) , COS(  X ) ) 
),  - S I N ( X ) ) 


ANR (N) 


i 

I 


r.l 


i 


\ 


f 


'WJRi  w»P!S"  ■ jpwh.w  ■"“" 


’TTTT’V'n’ 


ONE  = -1. 

no  to  n = i,  200 

ONE  = -ONE 
FN  = N 

Cl  = 2.  * FN  - 1. 

FN  = Cl  * EMI  / X - EM 2 

IF  ( N .LE.  NX  ) ANZ  = ANR(N) 

IF  ( h .GT.  NX  ) ANZ  = ANF ( N*  Z ) 

CFNOX  s FN  / x 

Ct  s HE*L(  t N ) 

C?  « WEAK  EMI  ) 


A !<l  3 

( 

( ANZ 

/ 0 ♦ 

CFNUX 

) 

* 

Cl  " C2 

) / 

1 

( 

( ANZ 

/ 0 + 

CFNOX 

) 

* 

FN  - EMI 

) 

BN  s 

C 

( 0 « 

ANZ  ♦ 

CFNOX 

) 

* 

Cl  - C 2 

> / 

1 

( 

( 0 * 

ANZ  ♦ 

CFNUX 

) 

* 

EN  - EMI 

) 

XF  ACT 

3 

z.  * 

F N ♦ 1 

• 

XFACB 

3 

ONE 

* ( F N 

T O.S 

) 

oe*t 

3 

UF.XT 

♦ XFACT 

* ME  AL  ( 

AN  T BN  ) 

OSC  AsosC  A + XF  ACT»  ( C ABS  C AN  )**2+CAtiS(RN  )**«?) 

0B3C  = UBSC  ♦ XFACH  * ( AN  - BN  ) 

EM?  s FM1 
fnl  s FN 

IF  f w .td.  1 .OH,  E N ,LT.  t.2  * X ) GO  to  ?0 
IF  ( AflS  ( t.  - AM  A X 1 ( NEXT,  OEXTST  ) / AMJNK  Ofc  X T , 
1 .It.  l.E-S  1 Oli  TO  <J0 

20  OEXTST  * NEXT 
30  CONTINUE 

00  QEXT  s 2.  * OF  X T / X2 
(JSCAs?.0*uSCA/X2 

OHS  x a.  * ( CAHSl  08SC  ) **  2 ) / X2 

HtTUHN 

EnO 


QtXTST  ) ) 


-i 

l 


l 

i 


t 


j 


1 


I 


SimKillI  r 1 Mt  Mjj.  ( X,  |,«,  II,  QSfA,  MxT,  s ) 

TmIs  IS  » Mlt  Pf  III  T I NF  - *"  lOE  t FijIJAT  lll'iS  Apt  *)L'C  UMF  N T F T N 

G E 77TMP-?,, 

This  WOuTlNt  uses  M{F  ThEuPV  Ti  I C a LCiJI.  ATE  THfc  E F t-  l C J E >iQ  I *■  S F . 

SC4TTFP1M,  AMn  aPS'IPPT  JON  AM)  ThF  St  A T T E P I M,  KATTFP*.  EOk  A 

ijn iFnpM  shhemcaL  PAFTTOLE 

TSPljT  s 

X = M1PMAL  T7Fn  SIZE  PAPA*«FTFk,  *"■  n I C P F 1 1I J A L S I*H,  Pi  T I -F  S 

The;  pauius  he  the.  spmfke  r- i v t ri f < 1 hy  the  havElF'iGTh  m im 
Incident  paDIatI'.iaj 

OP  - PE  A L PAPT  if  The  COMPLEX  I nDF  x i -F  hEFNaCIIon  ' f-f 

jPhFkE 

I>1  z IHAr-INAtiy  EA«7  f)F  THY  COMPIFx  [M'Fx  I’F  F E F ►>  At.  I l ' 1 ■ «F  Iff 

sPhf pf 

MiTE  THAT  ThP  CIME’LEx  I'jDEX  )>E  HE  F <<  A C T I ijl-  IS  aSS",,E|'  I : n 

N ; HP  - I * | » T 

nuTKi  I s 

'•JSCA  = SCATTEwIl.r;  FFEIC1F'-Cv,  »HlCM  F u 1 1 A l.  S 1 HE.  jCATTFrTm,  r»<ss 
SECT  J,  IN  IF  THE'  SPMFhf  LIvJoF.0  hY  The  C*-"S^  SF  f 1 I • "J ».  I jh  < 
I iF  ThF  SHhEPE  ( S I bMA  / ( p 1 *p  AD  T 1 1 S * *P  ) 

Qfcxl  = t x TTnCTTL'N  EFFICIENCY,  *hICh  FuI'aI.S  Tmf  T f ■ T a l r SC  a I H i-  I •• 

♦ A F<  fjnx-  P T I » <»•  ) CKUSS  SECTION  i iF  1 he  SPMi.f  0 1 V I r F I.  h,  r~E 

CPn.SS  SECTIONAL  APEA  f-E  THE  sPhEM 

S = SCu  TTFE  T NC,  P A T T E k A'  OF  THE  hAhTaTJoU  SCaTTFsEh  -Y  ThA 

SPhFPE  , ASSU'-'ING  JnCTOEnI  iiM’iilau  t/E"  r A 1.1 1 A r I'  .'I  , Sl.ll  s 
•SCaTTEPInT,  function  Fill-  THfc  SCaTTE.  H I vF  A-.r,LE  ►ih.'SF  C.sjmF 
IS  O.IMJ-M),  S TS  iiM-Ni  )P*-i  A L I ZE  o , T 1 1 a 1 1 S , fit  r'TEU.-ai 

OF  S Ovfc>  a •=!  STEP  An  Tans  E'JUAi.S  PI*'.SCA*x**p 


IT  J M E N S I LI  N Sf.Pll,  XHilfpi),  SP(rM),  PPIIPI),  PPpIPII, 

i pTicpn,  PT?c?n,  pTt?n,  pp<?n 

r> i «e n.sti in  anP(?oo) 

Complex  0,  7,  E^l,  o*,),  Fn,  s,,  s,-,  a.f 
Cl  I M P | F X A.,|P,  f_  v,  ANZ,  an,  hn 
Cf'iPLF*  Ca,  Cl,  Cc,  CFnox,  C a F * C T 

(jaTA  x mu  / -1.,  -,o,  -.*J,  -.7,  -,s,  -.S,  -.j,  -.S,  -,d,  -.1 
l # 1 , » , A , t S f ,0,  . 7 , . P , 1 ( / 


SET  VALUES  oF  p x Index  of  kEFSactJon,  SIZE  pfpi-MFp,  a? 

m j e VAPIAMLF 
0 s r M p L * ( ( P,  -11!  ) 

C x z x 
7 z C«  » o 

SF  T II  T T I A L V A L • I F S OF  » I C C A T I .HE S SF  l FUNCTION 
FM1  = CM°l*(  SI'vf  X ),  cost  x 1 ) 
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o n o o 


r 

( 


*■«? 


rr^LM  Cf'S(  X 1 


SIN(  X ) ) 


I 

! 


C 

c 

r 


ZfWM  .11(1  FFFJCIFNCV  VAPIAhI.F  S AND  Sr*TTf  M‘‘l’  KAlTxPN  VAh]f>'Lf 
III.  S I S 1 , ?! 

sim  s o. 

SPM1  S l:  . 

ppi  m = «. 
pt  i ( n = ii. 

PPP<  T ) 

P1P( t ) s 0. 

S Ci'MTTM't 
QSfA  : I). 
f!t*l  = >*. 

<;fT  "p  A O w A Y ijF  *’<f  )/  A | i jf  ,■>  - ' I SF  I.  F *<  T / f>  A f K ».  A ■«  f)  S - F C '.I I f1 

Tt  C"*<1  MU*- 

Nx  s 1 ,S  ♦ X 

NX  s »-A<0(  P,  M T N 0 ( NX,  Poo  1 1 
«rgb  ( > , « 1 s A t.F  ( '■<,  / ) 

N I •*  1 = NX  - 1 

nil  t n T = 1,  N x k 1 

N S f * * 1 - J 

C N = F L 1 1 A t ( N 1 

ANPfN-1)  = C*J  X 2 - (1,*0.)  X f CN  X / ♦ 1 } 

10  C ‘ IN  T t V It 


CAtC"|.ATt  TiF  FFFICIFNCTFS  A(r  ST  A T H *•  T » r>  °ATlFt->.  HaI*‘J  T“iF 
INF  I M Tf  nF_PTFS  FXPA'JSI'I*  FiiFMII  «s 

1)1  1 1 K 3 1 , p II  -I 

FI  8 N 

C 1 S i . * FI.  - 1 . 

FA.  s Cl  * F * 1 / r»  - F * P 

if  f , .It.  “I « 1 AN/  = AN*  (M 

IF  ( v , T,  T.  NX  ) AN/  = t»’l  f / ) 

CFNMX  s Fft,  / t 
Cl  = -f»t ( F \ 1 

CP  A PF  AL ( F M I ) 


AN  : 

( ( an? 

/ 

r*  4 

C F I lit 

) * r t - Cp  ) 

/ 

) 

( ( an/ 

/ 

t»  4 

r.  f N'fix 

U. 

• 

♦ 

l 

p-j  - 

( ( <)  * 

AN/  ♦ 

CFMix 

) » r i - f p ) 

/ 

1 

( ( U * 

A 

N/  ♦ 

C.  F Mi* 

) . F-  - FM 

) 

XF  Af. 

T = f U. 

* 

F *■.  ♦ 

?.  ) 

/ » ♦ »p 

NSC  A 

= i^c* 

4 

XF  AC  T 

» ( 

Caps*  a i j • » p 

♦ 

TF  X T 

8 |JF  x 1 

4 

*F  AC  1 

4 H t A L ( A \ 4 *sN  ) 

fc*P 

3 F"  1 

F-l 

s FN 

Ol*  u 

I s 1 » 

?\ 

IF  ( 

N . 1*  r , 

p 

1 r,,i 

in  t, 

IF  ( 

'■  .k'J. 

p 

1 r,._; 

T IF  7 

ppm  s i . 
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1 
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p r ( i ) = «mij(  1 1 

Go  to  » 

7 ppm  = 3.  * x«u( i ) 

pirn  = 6.  * x”u( n **?  - i. 

Gu  to  « 

ft  PP(I1  S f < 2.  * E N - 1 . ) * X«!J(  1 1 * PP  1(1)  - y K'  * PPPCJ)  ) / 

i < y m - i . ) 

prm  = x*ud)  * < pp  c i ) - pp?M)  ) - < ?.  * Fiv  - t . > * 

1 ( 1.  - xmU(I)  **2  ) * ppirn  ♦ p T?  f I ) 

8 CXFACT  s ( 2.  * FN  ♦ 1.  ) / ( FN  **2  ♦ FN  ) 

Cl  s PP(I) 

C 2 S P T ( I ) 

Sim  S SKI)  ♦ CXFACT  * ( AN  * Cl  ♦ BN  * C2  ) 

S2 ( I ) s S2 ( I ) ♦ CXFACT  * ( BN  * Cl  ♦ An  * C2  ) 

PP2 ( I ) = PP1 ( I) 

p t? ( n s prim 

PP 1 ( T ) S PP ( I ) 

PT1  (I)  s PTU) 

« CONTINUE 

I F ( n ,£Q.  1 .OP.  FN  ,LT.  1.2  * X ) GO  IU  2 
CHECK  IF  THF  INIFINITE  SERIES  has  CONVERGED 

IF  f A R S ( 1.  - AMAXK  QEXT,  OFXTST  ) / AMIN1C  GMT,  OEXTST  ) ) 

1 .LE.  S.F-3  ) GU  TO  3 

CONVERGENCE  HAS  NOT  BEEN  REACHED,  COMPUTE  NEXT  TERM  IN  SERIES 
2 USCAST  = USCA 
OEXTST  r NEXT 
1 CONTINUE 

SERIES  HAS  CONVERGED,  COMPUTE  THE  UNPOLARIZED  SCATTERING  PATTERN 

3009  I = 1 , 21 

S(I)  = 0.5  * C CABSI  S1(I)  ) **2  ♦ CABSC  S2(I)  ) **2  ) 

9 CONTINUF 


RE  TURN 


onn  n n n n on  no  nnnnnnnnnnnnnnnnnnnnnnnnnnnnn 


COMPLEX  FUNCTION  ANF ( INDEX*  Z ) 


iir. 


THIS  15  A WOE  ROUTINE  - DOCUMENTED  IN  REPORT  GE77TMP-22 
THIS  ROUTINF  IS  CALLED  «Y  MISCAT 

THIS  FUNCTION  EVALUATES  THE  COMPLEX  QUANTITY  A(N,Z)  WHICH  IS  USED 
IN  THE  MIE  FORMULAS*  WHFRE 
A ( N » Z ) = -N/Z  ♦ J(N-1/2*Z> /J(N+1/2*Z> 

Z = M*ALPHA 

N = M(PEAL)-I*M( IMAGINARY)  = COMPLEX  INDEX  OF  REFRACTION 

ALPHA  = 2*P I »R/WAVELENGT K = NORMALIZED  SIZE  PARAMETER 
R = RADIUS  OF  SDHEPE 

N = ORDER  OF  THE  FUNCTION 

J = BESSEL  FUNCTION  OF  COMPLFX  ARGUMENT  AND  HALF-INTEGER 

ORDEP 

THE  METHOD  OF  EVALUATION  USES  THE  CONTINUED  FRACTION  ALGORITHM  OF 
WILLIAM  J LENTZ  - GENERATING  BESSEL  FUNCTIONS  IN  MIE  SCATTERING 
CALCULATIONS  USING  CONTINUED  FRACTIONS 
APPLIED  OPTICS*  VOL.  15.  NO.  3*  MARCH  1976 


INPUTS 

INnrX  = ORDEP  OF  A ( N * Z ) * THAT  IS*  INDEX  = N 
Z = COMPLEX  ARGUMENT 

OUTPUT 

ANF  = A ( N * Z ) 


COMnLFX  Z*  N * D,  T,  PN,  PD,  Tl,  T2  , E 

DEFINF  ARITHMETIC  STATEMENT 

C<  X > = 2.  * S * ( FN  - 0.5  ♦ XI  ) 

SET  VALUE  OF  FIRST  PARTIAL  FRACTION  TERM  FOR  NUMERATOR  (PN) 

FN  = INDEX 
S = -1. 

CP  = 2.  * FN  ♦ I. 

PN=CP/Z 

SET  VALUE  OF  FIRST  PARTIAL  CONVERGENT  FOR  NUMERATOR  ( N ) 

N = P'I 

CALCULATE  SECOND  PARTIAL  FRACTION  AND  CONVERGENT  FOR  NUMERATOR 
CP  = -2.  * FN  - 3. 

T=C P/Z 

PN=T ♦ ( I . ,0, ) /PN 
N=N#PN 

SET  VALUE  OF  FIRST  PARTIAL  FRACTION  (PD)  AND  CONVERGENT  (D)  FOR 
DENOMINATOR 


101 


• W4  v 


PD  = T 
D = P 0 


CALCULATE  THE  HIGh£R  ORDERS  OF  THE  PARTIAL  FRACTIONS  AND 
CONVERGENT  S 

X I =2  • 

DO  TO  J = 1 , 100 

X I = X T ♦ 1 . 

s = -s 

T = Ct  X ) / 2 
PN=T* (1. ,0. ) /PN 
PD=T+(1.,0.)/PD 


IN  THf  RARE  INSTANCE  THAT  THE  NUMERATOR  PARTIAL  FRACTION  TERM  IS 
NEAR  ZERO « USE  THE  LENTZ  ALGORITHM  IMPROVEMENT  METHOD  TO  INSURE 
ACCURACY 

IF  ( CABS ( PN  ) .GT.  l.F-A  > GO  TO  20 
S = -S 

XI  = XI  ♦ 1. 

T 1 = C ( X ) / Z 
£ = T 1 * PN  ♦ ( 1 . ,0.  ) 

N s N # E 
S = -S 

XI  = XI  ♦ 1. 

T2  = Ct  X ) / Z 
PN  = T2  ♦ PN  / F. 

IF  THE  DENOMINATOR  PARTIAL  FRACTION  TERM  IS  NEAR  ZERO.  USE  THF 
ALGORITHM  IMPROVFNENT  METHOD 
IF  ( CABS ( PD  ) .GT.  l.F-A  ) GO  TO  10 
E = T 1 * PD  ♦ ( 1.  .0. ) 

D = D # E 

PD  = T2  + PD  / C 

GO  TO  20 

10  D = D * PD 

PD=TU(1.,0.)  /PI) 

D = 0 * PD 
PD=T2* ( 1 . .0. ) /PD 


20  N = N * PN 
0=0*  PD 

CHECK  IF  CONVERGENCE  HAS  BEEN  REACHED 

IF  ( ABSt  CABS ( PN  ) / CABS!  PD  ) - 1 . ) .LE.  l.E-6  ) GO  TO  AO 
30  CONTINUE 


CONVEPGENCE  HAS  BEEN  REACHED,  SET  VALUE  OF  ANF 
AO  ANF  = -FN  / Z ♦ N / 0 

RETURN 

END 


I 


APPENDIX  B 


IMPLEMENTATION  OF  THE  GENERALIZED  POWER  LAW  SIZE  DISTRIBUTION 


The  power  law  size  distribution  for  nuclear  produced  dust  particles  ; 

is  given  by  Equation  26  t where  p is  the  power  law  exponent.  In  the 
present  WEPH  model,  p is  taken  as  4.0,  a value  representative  of  dust 
particles  generated  from  loose  unconsolidated  soils.  The  present  fixed 
p model  is  easily  generalized  to  an  arbitrary  p.  We  first  present  those 
model  equations  which  are  changed  due  to  an  arbitrary  exponent;  then  we 
give  a computer  listing  of  the  revised  computer  routine  PGROUP.  PGROUP  l 

calculates  the  extinction  and  backscatter  cross  sections  per  particle 
for  a given  size  interval  of  dust  particles. 


The  generalized  model  equations  are: 
Number  distribution 


f (a)  = Ks  a'P 


Total  number  of  particles 
k 


N 


PT  p - 1 


1-p  1-p 

a - a„ 
s x. 


Fraction  of  particles  in  size  group  i 


P X 4 
P = 4 


(B-l) 


i 


(B-2) 


(B-5)  I I 

I I 


(B-4) 
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w 

Ci 


Normalization  constant 


K = 
s 


5.44  x 1012  R N 



"P,. 


Backscatter  cross  section 


J 

4 (a1  ~ ^ 
i i+l;  j=l 


a - 10  rc(P-P 

bi 


4 - |> 


4-P  4-p 

a - a 

s 


a£ 
£n  — 
a 

s 


bj 


P ^ 4 


P = 4 


(B-5) 


(B-6) 


exp  (xb^+3-p)«.n 


a . 

i_ 


o,  . = K,  (a . ) 
bj  b^ 


p-o,  , , 

aj  ^bj^ 


a • 


— — An 
x.  . a . 

i.bj  J 

J 


l.xtinction  cross  section 


-4  J 

10  7T  (p  - 1 J 


O 1 


0 nP-lJ  V rr 

1-i'  .l-p,  Lj  °ej 


4 (a  1 - a1  ) “ 
i i+lJ  j=l 


exp  I (x^+3-P)  £n 


^1-1 

3j  j 


xbj  * <’  -5 


bj 


a.  (x  . +3- 11) 
J ej  ’ 


XCJ  * I’-3 


. K !.i  ) 
l’ i e y 


in 


a . 
J 


ej 


x . = p-3 
eJ 


,e  computer  iistinj;  of  subroutine  PGROIIP  follows. 


(B-7) 


(B-8  J 


(B-9) 
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<5 1 1 « k r, , I T I Ni t pr.HdnP(  JGM'uP,  E~  hk n , 


I)  1 1 C n , U J t C I 


) 


C 

c 

c 

( 

c 

c 

c 

c 

r 

c 

c 

c 

r 

c 

c 

r 

c 

c 

c 

c 

c 

c 

r 

r 


c 

c 

r 


PAUL  POUTINF  him:  If  r f-  0 flip  It  S f Tn  wf^H 

this  H'HiriNf  computes  the  a vFpage  bace.sl a r re p ah> 


F X T i mC  T i Cif. 


O'I’SS  SECTIONS  FOP  A SPECIFIED  DUST  PAPIICLF 


INPUTS  FPOM  r&l.l  STaTEmENT 
FkflO  s FPtrailFNCV,  HH7 

Jf.kl'ip  s UUMHER  nf  DUST  PARTICLE  GPuOP 

niflk  s HEAL  PapT  (IF  r.lJi''PLE*  INDEX  OF  pFEPACTILn 

!>)cfT  = IMAGINARY  PANT  OF  CO-PLEX  I^uE*  IF  P t F P A (.  I J <n.  I ncTE  I m M 

THE  INDEX  OF  PFFPACTI'IU  IS  * s FIFO  - I * 0 1 F C I Si  i"Al 

both  iufcp  a no  oifl!  ape  Positive  ) 

P s tXPliNfriT  HE  G E N F p A L I 7 f 1 ' Eil^K  LAf  P Pi .P  A H I L J I Y 

DISTNIM'T  Ihn 


OUTPUTS  T r i TEST  APea 

Slot  = a / e page  extinction  cposs  section,  *v 
sigh  = a v E page,  paces c aTTEp  cposs  section,  »d 
s i g s = a v i page  scattep  cpuss  SFi.Tin.v,  v 

c i /Tt  ST/  NPG,  AKAP  (11  ) , SIGH  ( 1 ()  ) , SIGt  M (I  ) , S I G S U 0 ) 

INTTITALTZE  VAPIAhLES  Fijp  COMPUTING  avEpAGE  HACeSC.ATTEFv  aiD 
extinction  I.POSS  SECTIONS 

P r :I. I 4 1 SD,?7 
x I.  A M n A s i . E 4 / F P F 0 
ExsE’t/XLAHOA 
I PG  = jr,p.j'JP 

I pg  t = r °g ♦ i 

A1 =APAO(IPGI 

X = A 1 *E  X 

CALL  "ISCATC  x,  DIFCPx  DlECI*  xk«I,  xkE  T . xksi  ) 

AJsIHi-dEGI  1 
A 3s  A?** ( 1 . »P ) 

' 3 I sA 1 * * ( 1 . u-P  ) 

E NF 1 = ( A 3 I -A  3 1 
DEL  T AsAP-A 1 

jHAx=MAxn(l,S*IEIX(AHlNIM.0,(AP/Atl/jOT) 

DEL  T AsOF.L  T A /ELOA  T ( JHAX  I 
S G 4 = 0 , 

SGEsCT  . 

SgSsO  . 0 

Do  toy  J = 1 , .1 M A X 
AP  = A 1 + DEL  I A 

XsA^.E  X 

Call  “I  SC  ATT  x,  ojecp,  0 1 E C I , Xep,  xkE,  xcS  ) 

COMPUTE  P,j«f  U L4*  APPPOXIMATInN  E(|p  HACeSlATTEP  and  E X T 1 n C T I i ' » 
E E E I c I e uCt  S 
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c 


c 

c 


c 

c 


c 

c 

r 


c 


c 

c 

c 

r 


AL°f’  A = * LCk,  ( a?/ A 1 1 

K^zAt  or,  C XKH/ XKH  n / ALflliA  j 

xe z AL'ir, t x*^  / xk »- n / ai  itu« 

xS  = al  IGIxks/xkSI  )/al<ii,a 

i 

niMpiiTF  scattfp  cross  stCTiuK 

x 1 z x s-p» 5.0  1 

IF(  XI  ,FU'.  0.0  ) Gl'  Hi  <!  ] 

SuSJz«KSI*rPxP(xt*ALiir:A)-i.())/fx|*Ai**(P-5.o)) 

RO  TO  4 

? Sus.jzxKSl  *A|_nr,A/A  i **vs 
4 L • t K ' T T \ . II 

CuMPHTf  hACK^CATTtP  C Pi  USS  5M-TI0N 
X 1 z x«-P+  5 . 

T F I X n HI » 1 0 , 1 0 

in  <sGRJ=xKKT*(t*Pfxl  * AL  OUA  )-i.)/(xi*ai**(p-5.1) 

Gii  TO  4 0 

«»0  <?C-H  J=  x*h  I * ALOGA/ A 1 **y  H 

CuMPUTF  F x T 1 FiC  T Iiin  CPnSS  SECTION 
40  XlzxF-P+5. 

JPFxi  )So,e>o,S<>  .! 

SO  SgF Jzxkg 1 * (F xp ( » 1 * ALPbA )- 1 . ) / r X 1 ♦ a 1 * * (P-5.  ) ) i 

Gl:  m All  i 

60SGFJ=X*t1*Au)GA/Al»*xfc-  j 

ACC  I”|.ILA  TF.  CROSS  SFCTIuSS 

00  SSP  = <?r,N  + eGP.I  i 

SGFzSGF+SbFJ  1 

SGSsSl.S  + SbSJ  1 

A 1 z A P j 

XK  0 J = XK  pi 
X»F  I z X x t 
X k S J :*«S 
1 00  CUM  TNiifc 

Ci.iMPi.lfF  A y F P A G F ^ACK SCATTFP  »,ri  FxT  IM  TJiin  (.PUSS  SFf.TIuNS  Cm.sS 
SFCTI-'iPS  FFp  PAPTICLF  (m**?) 

SI  GO  f I PG 1 = 7 ,^S4E -b* ( P- t . 0 ) *SG«/F NP1 

STGtfIPbl  = 7.8S4t-S*CP-1.0)*SGF/F*JPl  J 

SlGSfl  Pb)s7.  RS4fc-S*(  P-1 .0  ) A.SGR/F  SP1 
PF  turn 

Fnc  : 


1 ....  • r i i -i  iti 
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APPENDIX  C 

SOLUTION  TO  CUBIC  MIXING  RULE  EQUATIONS 


The  two  mixing  rule  equations  are  (see  Section  4 ) 


3 e*e*2  + (3  e*2  - C)e* 


(C-l) 


e*-3  - 3 e*  e*2  + (3  e*2  - D)  e*  - e*3  = 0 
ni  dm  d md 


(C-2) 


where 


C = 


(C-3) 


D = 


(e*  - £*)• 

(l-4')3e* 

m 


(C-4) 


The  cubic  equations  (C-l  and  C-2)  can  be  solved  for  £*  and  £*,  respectively 
(the  other  parameters  assumed  known)  by  the  standard  cubic  equation  formulas. 
Some  care  must  be  taken  because  the  coefficients  are  complex.  Let 


i £*  for  Equation  C-l 

) 

\ e*  for  Equation  C-2 

(C-5) 

-3  £* 

J £d 

(C-6) 

2 

j 3 e*  - C for  Equation  C-l 

( 3 e*2  - D for  Equation  C-2 

(C-7) 

*3 

- C* 

£d 

(C-8) 

Our  cubic  equations  are  now  in  the  standard  form 
z3  + Pz2  + Qz  + R = 0 


(C-9) 


but  with  all  quantities  complex.  As  in  the  case  of  real  coefficients 
define 


■ i(3Q  • p2) 

* ^(2p3-91 


9PQ + 27R 


J b2  a5 

6 = VT+27 


f = - j b + e 


g = - y b - e 


A = f.  , the  first  complex  root  of  (f)' 


>83  = three  complex  roots  of  (g)' 


CC-10) 


(C-11) 


(C- 12) 


(C  - 1 3 ) 


(C - 14) 


CC- 15) 


(C- 16) 


it  does  not  matter  which  of  the  two  complex  roots  is  taken  for  e.  By  con- 
vention we  take  the  first  root.  The  three  complex  roots  of  f and  g are 
ordered  in  the  following  manner.  Write 


f = fR  + ifx  = rex"  , 


(C- 17) 


where 


J c2  c2 

r ■ % * fI 

i = /T 


-TT  < 0 < 7T 


CC  - 18) 


fD,fT,r,  and  6 are  all  real  quantities. 
K I 


The  three  roots  of  (f)  ' are  ordered  as 


, 1/3  i0/3 

f^  = r e 


(C- 19) 


f = rl/3  gi  (6/3  + 2/3ti) 


f. 


rl/3  ei(0/3  + 4/3TT) 


, 1/3 


(C- 20) 
(C- 2 1 ) 


The  three  roots  of  (g)  are  similarly  ordered.  Form  the  three  trial 
solutions , 


ll  “ A + *1 


t7  = A + g. 


t,  = A + g, 


(C - 22 ) 
(C  - 2 3 ) 
(C-24) 


Substitute  the  trial  solutions  into  the  reduced  complex  cubic  equation 


t + at  + b = 0 


(C  - 25 ) 


One  of  the  three  trial  solutions  will  satisfy  the  reduced  equation.  Let 
B be  that  value  of  g^,  g^,  g,  of  the  successful  trial  solution.  Then  th( 
three  complex  solutions  to  the  complex  cubic  liquation  C-9  are 


‘l  - A + B - i P 


- J (A  ♦ B)  ♦ (A  - B)  - | P 

- j (A  ♦ B)  - (A  - B)  - ip 


(C-26) 

(C- 27) 

(C-28) 


A computer  program  has  been  written  to  solve  Equations  C-l  and  C-2. 
The  inputs  for  Equation  C-l  are  e*  , and  tp;  for  Equation  C-2  the  in- 
puts are  e*  , and  (J).  The  input  format  is  1 1 ,E9. 0,4E10. 0 , and  the  in- 
put parameters  are  arranged  on  the  data  cards  as  shown  below: 


i 


Column  1 10  20  30  40  50 

1 Ue(c*j  I Im (e* } I Re(t*)  I Im(e;)  I <|>  I 


2 Rett*)  Im(e*)  Re(c*)  Im(e*)  4- 


Re  = real  part  of 


1m  _ imaginary  part  of 


The  1 in  column  1 indicates  that  it  is  Equation  C-l  that  is  to  be  solved; 
a 2 in  column  1 indicates  that  Equation  C-2  is  to  be  solved.  As  many 
cases  may  be  stacked  as  desired;  a blank  data  card  is  placed  at  the  end 
of  the  data  deck  to  signal  end  of  data.  Note  that  the  input  for  the 
imaginary  parts  of  the  relative  permittivities  are  negative.  That  is. 


Re(e*d)  = e’  > 0 

lm(e*d)  = - cj’  I 0 . 

The  output  consists  of 

1 • The  input  parameters 

2.  The  real  and  imaginary  parts  of  P,  Q,  R of  the  cubic  Equation  C-9 

3.  The  real  and  imaginary  parts  of  the  complex  solutions  Zl,  Z2,  Z5  of 
Equations  C-26,  C-27,  and  C-28 

4.  The  real  and  imaginary  parts  of  the  square  roots  of  Zl,  Z2,  Z3. 

One  of  the  three  solutions  Z^,  Z9,  Z3  is  the  proper  solution  for  z* 
for  Equation  C-l,  or  e*  for  Equation  C-2.  We  have  arranged  the  program 
so  that  for  all  cases  we  have  run,  the  proper  solution  has  been  Z^.  It 
is  conceivable  that  on  a different  computer  or  for  some  special  input 
values, the  proper  root  may  not  be  Zj.  Therefore  as  a precaution,  we  print 
out  all  three  roots.  Generally  the  proper  solution  is  obvious.  The 
proper  solution  for  e*  or  must  have  a positive  real  part  (c')  and  a 
negative  imaginary  part  (-e").  Two  of  the  solutions  may  satisfy  the 


r 


positive  real  and  negative  imaginary  criteria,  but  in  our  experience 
the  proper  solution  is  obvious  from  the  magnitudes.  The  square  root  of 
the  relative  permittivity  is  also  computed  since  this  is  the  index  of 
refraction,  which  is  needed  for  propagation  calculations, 


Note  that  m is  normally  written 


m = 


mR  ‘ lml 


(C-29) 


(C-30) 


so  that 


= Re  (m) 


(C-31) 


m.  = -lm(m) 


(C-32) 


As  in  the  case  of  c*  and  £*,  the  output  for  m has  a positive  real  part 
and  a negative  imaginary  part. 

The  listing  of  the  computer  program  follows.  The  program  consists 
of  a driver  routine  and  the  two  routines  which  solve  the  complex  cubic 
liquation  C-9.  The  complex  cubic  equation  routines  are  general  routines 
which  can  be  separated  from  the  mixing  rule  driver  and  used  to  solve  any 
complex  cubic  equation. 


&1 


E 


f. 

f. 


t 

f:  ! 


f 

f 

K 


i 


c 

c 

c 

c 

c 

c 

c 

c 

r 

C 

C 

C 

C 

r 

c 

c 

c 

r 

c 

c 

c 

r 

C 

c 

r 

r 

c 


c 

r 

c 


c 

c 

r 

c 


Tw,s  IS  THfc  fiPIVEP  PiRiTIME  FI'S  Tmi-.  COhPLEx  CUMC  "Ixlfr,  si.ik 

0')C"MFNT  AT  ION  is  GIVEN  fN  APPfM/1*  C 'if 

iHjst  ci  nun  poofling  am.-  pbupac.at  jom  effects 
Fits  k«r;»h  AMU  CHM^UMICAT  JhmS  Cofif  S 
UF7SfMp-81  OCT.lPtP  tP7e 

t’Y  JA'IF  S H THOMPSON,  GENERAL  FIFCTPIC  - (E-Pi; 


C'|Mpl  t » E » FM,  FT),  C,  P,  P,  <g,  P,  1 1,  7i 


IT4PF  Is  THF  [MPIJI  tape.  MijMhEP,  JTApF  IS  THE  nu|Pt;T  TaPF 
The  Values  have  hpe-,  set  at  ftape  = s a *0  JTaFf  = h m,i!  t.-f^f 

rtF  CHANGED  Flip  A ‘g  Y COHPuTfP  My  “F  SF  IT  I MG  The  F 1 1 (_  I i : a I N G ITAPt  * 
JT  APF  S r A TFMF  nTS 
ITAPF  = S 
J T APF  = t, 


V. 

r a . 


PE  A 0 TMpiil  FOP  F J p s T CASE 

10  PEAOf  T r A PE  , 11  1 POPE,  Vp),  *n,  xp?,  x T ? , PM  I 

11  fubmatf  II,  F<g.O,  UtlU.G  ) 

Check  TF  all  CASES  have  hefn  bum 

IFl  PQDE  ,LF. . 0 1 STOP 

SET  VA(  II?  OF  E0 

Fn  = c hpl * ( VP i , x 1 1 ) 

Check  *mTCH  Efgl'ATIflN  Is  Tl,  riE  SOLVED,  C-1  OP  c-i 
IF(  -line.  .Fd.  ? 1 go  rn  ? o 


E 0 1 1 A T I Cl  N c-1  IS  TtJ  HF  SOLVED 
E«  s CMPL  A ( VP?,  *1?  ) 

C s ( 1.0  - PmT  )**1  * ( EH  - F 0 )**?  / Fh 

n s s.  * tn**2  - c 

Go  TO  so 


Equation  c-?  IS  Tf)  HE  SuLVFO 
?0  E = CMPL  V ( XP?,  xl?  1 

0 : ( F . ED  1**3  / ( ( 1.  - PHI  )*»}  • e ) 

Q = 3.  * tO*  *?  - [) 


10  P s -i.  a ED 
P r -E  0 * * i 

CALL  CCUHICf  p , (J , P , 7 1 , 2 ? , ^ 1 ) 


1 


i 


r 

r 


"MTf  f M t T [iPiiT  , P-i.-p  VtLHcS,  A ' ■ 1 1 T *->  H IHUj-f.  (^"'Pl-FA  >-  < t • T S > i 1 * F 
M I * I M,  PHI  f F I*  T 1 O '. 

If  ( “HOF  .Mj.  ? i r.f  r ' i Sn 


I 


r 


b 


r 


c 


r 


C 


r 


r 

r 


c 


c 

r 

f 


"&ITFI  JTapf,  up  ) F-nriF  , *(»|,  » I i , k]P,  Pf-i 

»0  F i iO"A  T { 1 mu,  M "F  u Im'Fi  ofAl(FO)  J ■«  A M r a " y ( F u ) 

I I " * I,  T \ A 0 Y ( F « > v..ILIimF  fw»CII'l*<  / 

P I " , 1 S , 5 * , 1 PSr  1 4.  n 

fin  T ■ i hi) 

50  *01  TF  ( JTAPF  , 5 1 ) "nuF  , <H,  YJ1>  *>*?,  KIP,  l-il 

51  fiip-a  r f i ho  , h i hf  o i Muf  > -fal(fo)  i»  a i-  t f a«y  (F  •• } 

1 I 4 A (;  T M A P Y ( F 1 Vllll'Pf  FPACIT"'  / 

PI  " is,  it  . 1PSF  1 4.  SI 
CO  *o  l , F ( JTapf  , el  ) 

ol  FilP"Aff1*"4,SSP 

IMS  / 

? 1 " , 7 P 4 0 1 A L ( P 1 

5 I 'A  A r,  I A.  A « Y f K > / 

4 1 h , t u*,f  i 4 . ;s  i 

IF  ( M |fF  . F U . P 1 (ill  T f • M> 

*“llf(  JTaPf,  7,1  ) /I,  /P,  73 

70  F , iO~  AT  ( | p'j,  7f>H  TmpFF  POL  rs  ,.'F  CIPPLF*  C'-"IC 

1 M I * I OUlfc  F (HATlUf.  f-t  / 

PI"  ,»S"  (11  Pf  A ( ( F 1 1 "A',  If, AS  Y (F  1 (?>  »FM.(F)  (pV(F  ) 

5 ( 5 ) PF  AL  ( ► 1 TsAIiInAmV  (F  ) / 

PI P .1 PSfc 1 4. J 1 

r,(,  ir.  ao 

KO  p p T I F I JTAPF , PI  ) 7 1 , 7?,  73 

" 1 F I iO  v A T t 1 h,,,  7f>w  T ><»F  F kHi.TF  ,F  CifiFI.Fx  f.  i . • 1 C 

1“lM‘f.  PUt  F F pu  a T l . H r-P  / 

P 1 U , P S P ( 1 ) PF  A L ( F M ) I P A G 1 V A t-  Y ( F P I I 1 PF  A I f F ‘ ) I >'  A (.  1 !.  f.  fi  Y ( r • ) 

5 ( J 1 PtAu(F“)  ]"AGINaOY(F"1  / 

41"  , IPSt 14. 51 

C,)MPnrF  ISUFX  <lF  PFfPArilfiS 
40  71=15^01(711 

7p  = c 5 JP T ( 7 P 1 
73s(.5  WT  f 73) 

"01  IF  ( JTaPF  , Q 1 ) 71  , 7?.  I 5 

4 1 Fi  Pp  1 A r f 1 HO  , S 1 rl 

P 1 F-  ,»S"  (1)  PFA(_l")  I " A I,  1 v A P < ( P ) 

5 (5)  ofaL(m)  r«lGlMPY("|  / 

51"  ,1F*>F.  14.5////1 


I VI F » 'IF  p F F p A t 1 1 I / 
(P)  OfcAL  ( ")  |p«i,I'*»M»  ) 


p > f , 4 

SI  ANl  A F 0 CUHir  F l.-UAT  luf.  tl'FFFIClr 
I“AG[AAF»(p)  w t A | f , j > IMAGINAPl  (f)  P F “ L t h ) 


P F A L.  I t u ) 


P F’  A I ( F ) 


fin  T n in 


i 


F if 
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SUBROUTINE  CCUBICt  P.0,S,Z1.Z2*Z3  ) 


f- 

jv 


THIS  ROUTINE  CALCULATES  THE  THREE  ROOTS  OF  A CUBIC  EQUATION  WITH 
COMPLEX  coefficients 

THE  CUBIC  EQUATION  IS  ASSUMED  TO  BE  OF  THE  FORM 

Z*#3*P»Z*»2+Q*Z*S=0 


INPUTS 

P = CONSTANT  COMPLEX  COEFFICIENT  MULTIPLYING  Z**2 

o = constant  complex  coefficient  multiplying  z 

S = CONSTANT  COMPLEX  COEFFICIENT 
OUTPUTS 

71  = FIRST  COMPLEX  ROOT 
Z2  = SECOND  COMPLEX  ROOT 
Z3  = THIRD  COMPLEX  ROOT 


COMPLEX  P.Q«S>Z1.Z2.23,A,B,C«D»E.F,T 

c 

A= ( 3.0*Q-P*P) /3.0 
B = (2.0*P*P*P-9.0*P*Q+27.0*S) /27.0 
F=CSOPT (B*B/4.0+A*A*A/27.0) 

F=-B/?«0+E 
X=REAL (F) 

Y=AIMAG(F) 

CALL  POLAR ( X « Y *R  » THET  A ) 

R=R**(l./3.) 

THETA=THETA/3.0 

CsCMPLXI  R*COS(THETA)  , P#SIN(THETA)  ) 

F = -B/ 2 • O-F 
X=REAL (F) 

Y=AIMAG(F) 

CALL  POLARt  XiY 1P1 THFT A ) 

R = R*»  ( I » /3. ) 

THFTA=THETA/3.0 

D = CM°LX  t P*COS(THCTA)  , R*S I N ( THFT  A ) ) 

Z I =C*D 

T = Zl*Zl*7.  WA*Z1  + B 

I F ( CABS(T)  ,LT.  l.E-3*AMAXl(  CABS(A*Z1>  . CABS(B)  ) ) GO  TO  10 
PI23=8.0*ATAN ( 1.0) /3.0 

DsCMPLXl  R*COS(THETA  + PI23)  . R*SI N ( THET A + P I 2 3 ) ) 

Z1=C*D 

T=Z1*Z1*Z1*A*Z1*R 

I F ( CABS(T)  .LT.  I . E-3* AMAX I ( CABS<A*Z1>  * CABS(B)  ) ) GO  TO  10 
PIA3=2.0*PI23 

I.'  = CMPLX(  P*C0S(THETA*PIA3)  » R*S!NtTHETA  + PlA3)  ) 

Z1=C*D 

T = ZI»Z1*ZUA»Z1*B 

I F t CABS(T)  .LT.  l.E-3*AMAXl(  CARS(A*Z1)  « CADS(B)  ) ) GO  TO  10 


SUBROUTINE  POLAR ( X«Y»R*THETA  ) 


T 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  ROUTINE  CONVERTS  a COMPLEX  number  FROM  RECTANGULAR (Z=X+IY) 
TO  °OLAR (Z=R*£XP(I*THETA)  ) NOTATION 


INPUTS 

X = REAL  PART  OF  COMPLEX  NU^BFR 

Y = IMAGINARY  PART  OF  COMPLEX  NUMBER 

OUTPUTS 

R = ABSOLUTE  MAGmjtudE  OF  COMPLEX  NUMBER 

THETA  = POLAR  ANGLE  OF  COMPLEX  NUMBER (RADIANS) 
THETA  IS  BETWEEN  -PI  AND  ♦ p i 


R=SORT(  X*#2  * y*#2  ) 

IF(  X , EQ.  0.0  ) GO  TO  10 

IF<  ABS(Y)  ,GT.  ARS(l.E12*X)  ) GO  To  10 

THETA=ATAN ( Y/X ) 

IF(  X ,GE.  0.0  ) RETURN 

THFTA  = TI-ETA*SIGN(1.0,Y)*<,.0«ATAN(1.0) 

RET'/RN 

10  The T As  2 «0» A TAN (1.0) 

IB(  V .LT.  O.C  ) TuETAs-THETA 
RETURN 

END 
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